{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "2.SimpleODE.ipynb",
      "provenance": [],
      "collapsed_sections": [
        "Y61YA90-WIZ1",
        "JksBmhyho5uZ"
      ],
      "machine_shape": "hm",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/jdtoscano94/Learning-PINNs-in-Pytorch-Physics-Informed-Machine-Learning/blob/main/2_SimpleODE.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Libraries"
      ],
      "metadata": {
        "id": "Y61YA90-WIZ1"
      }
    },
    {
      "cell_type": "code",
      "source": [
        " ! pip install pyDOE"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WZrPCr9GWEAJ",
        "outputId": "652e2d60-044a-4732-e41b-3e6f9d4e9a43"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
            "Collecting pyDOE\n",
            "  Downloading pyDOE-0.3.8.zip (22 kB)\n",
            "Requirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (from pyDOE) (1.21.6)\n",
            "Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (from pyDOE) (1.4.1)\n",
            "Building wheels for collected packages: pyDOE\n",
            "  Building wheel for pyDOE (setup.py) ... \u001b[?25l\u001b[?25hdone\n",
            "  Created wheel for pyDOE: filename=pyDOE-0.3.8-py3-none-any.whl size=18184 sha256=e5ea35a00c121de49f2b7a578cf68eedf4a8a2923a1bda18dbf1a9aec401cc9c\n",
            "  Stored in directory: /root/.cache/pip/wheels/83/ce/8a/87b25c685bfeca1872d13b8dc101e087a9c6e3fb5ebb47022a\n",
            "Successfully built pyDOE\n",
            "Installing collected packages: pyDOE\n",
            "Successfully installed pyDOE-0.3.8\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import torch\n",
        "import torch.autograd as autograd         # computation graph\n",
        "from torch import Tensor                  # tensor node in the computation graph\n",
        "import torch.nn as nn                     # neural networks\n",
        "import torch.optim as optim               # optimizers e.g. gradient descent, ADAM, etc.\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.gridspec as gridspec\n",
        "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "import matplotlib.ticker\n",
        "from sklearn.model_selection import train_test_split\n",
        "\n",
        "import numpy as np\n",
        "import time\n",
        "from pyDOE import lhs         #Latin Hypercube Sampling\n",
        "import scipy.io\n",
        "\n",
        "#Set default dtype to float32\n",
        "torch.set_default_dtype(torch.float)\n",
        "\n",
        "#PyTorch random number generator\n",
        "torch.manual_seed(1234)\n",
        "\n",
        "# Random number generators in other libraries\n",
        "np.random.seed(1234)\n",
        "\n",
        "# Device configuration\n",
        "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
        "\n",
        "print(device)\n",
        "\n",
        "if device == 'cuda': \n",
        "    print(torch.cuda.get_device_name()) \n"
      ],
      "metadata": {
        "id": "OdYMzuSDi7Js",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "a4871ba7-cbb0-485f-84d5-85a84b49c23b"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "cuda\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Tunning Parameters"
      ],
      "metadata": {
        "id": "JksBmhyho5uZ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "steps=5000\n",
        "lr=1e-3\n",
        "layers = np.array([1,50,50,20,50,50,1]) #5 hidden layers\n",
        "min=0\n",
        "max=2*np.pi\n",
        "total_points=500\n",
        "#Nu: Number of training points (2 as we onlt have 2 boundaries), # Nf: Number of collocation points (Evaluate PDE)\n",
        "Nu=2\n",
        "Nf=250"
      ],
      "metadata": {
        "id": "gL-vSDVLoyV7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Problem Setup\n",
        "\n",
        "### **ODE**\n",
        "\n",
        "$$y_x=\\frac{dy}{dx}=cos(x)$$\n",
        "\n",
        "The residual will be:\n",
        "\n",
        "$$f=y_x-cos(x)$$\n",
        "\n",
        "This function will be minimized by the NN\n",
        "\n",
        "**Initial Conditions:**\n",
        "\n",
        "$$y(0)=0$$\n",
        "\n",
        "$$y(2\\pi)=0$$\n",
        "\n",
        "\n",
        "**Solution:**\n",
        "\n",
        "$$\\int dy=\\int cos(x)dx$$\n",
        "\n",
        "$$y=sin(x)+C$$\n",
        "\n",
        "Since: \n",
        "\n",
        "$$y(0)=0=sin(0)+C$$\n",
        "\n",
        "$$C=0$$\n",
        "\n",
        "So: \n",
        "\n",
        "$$y(x)=sin(x)$$"
      ],
      "metadata": {
        "id": "koofMp-5Y-UA"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Functions"
      ],
      "metadata": {
        "id": "K40IpAXPjS4I"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def f_BC(x): # This function satisfies the boundry conditions. The same as the real one (To ease the data generation), but we may not have it.\n",
        "  return torch.sin(x)\n",
        "def PDE(x): # The PDE equation. We use it to get the residual in the Neurl Network.\n",
        "  return torch.cos(x)\n"
      ],
      "metadata": {
        "id": "hZFwvWyGs2RE"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Neural Network"
      ],
      "metadata": {
        "id": "ZHYqtyYJs3Jv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "class FCN(nn.Module):\n",
        "    ##Neural Network\n",
        "    def __init__(self,layers):\n",
        "        super().__init__() #call __init__ from parent class \n",
        "        'activation function'\n",
        "        self.activation = nn.Tanh()\n",
        "        'loss function'\n",
        "        self.loss_function = nn.MSELoss(reduction ='mean')\n",
        "        'Initialise neural network as a list using nn.Modulelist'  \n",
        "        self.linears = nn.ModuleList([nn.Linear(layers[i], layers[i+1]) for i in range(len(layers)-1)]) \n",
        "        self.iter = 0\n",
        "        'Xavier Normal Initialization'\n",
        "        # std = gain * sqrt(2/(input_dim+output_dim))\n",
        "        for i in range(len(layers)-1):\n",
        "            \n",
        "            # weights from a normal distribution with \n",
        "            # Recommended gain value for tanh = 5/3?\n",
        "            nn.init.xavier_normal_(self.linears[i].weight.data, gain=1.0)\n",
        "            \n",
        "            # set biases to zero\n",
        "            nn.init.zeros_(self.linears[i].bias.data)   \n",
        "    'foward pass'\n",
        "    def forward(self,x):\n",
        "        if torch.is_tensor(x) != True:         \n",
        "            x = torch.from_numpy(x)                \n",
        "        a = x.float()\n",
        "        for i in range(len(layers)-2):  \n",
        "            z = self.linears[i](a)              \n",
        "            a = self.activation(z)    \n",
        "        a = self.linears[-1](a)\n",
        "        return a\n",
        "    'Loss Functions'\n",
        "    #Loss BC\n",
        "    def lossBC(self,x_BC):\n",
        "      loss_BC=self.loss_function(self.forward(x_BC),f_BC(x_BC))\n",
        "      return loss_BC\n",
        "    #Loss PDE\n",
        "    def lossPDE(self,x_PDE):\n",
        "      g=x_PDE.clone()\n",
        "      g.requires_grad=True #Enable differentiation\n",
        "      f=self.forward(g)\n",
        "      f_x=autograd.grad(f,g,torch.ones([x_PDE.shape[0],1]).to(device),retain_graph=True, create_graph=True)[0]\n",
        "      loss_PDE=self.loss_function(f_x,PDE(g))\n",
        "      return loss_PDE\n",
        "      \n",
        "    def loss(self,x_BC,x_PDE):\n",
        "      loss_bc=self.lossBC(x_BC)\n",
        "      loss_pde=self.lossPDE(x_PDE)\n",
        "      return loss_bc+loss_pde\n",
        "\n",
        "    def closure(self):\n",
        "        \n",
        "        optimizer.zero_grad()\n",
        "        \n",
        "        loss = self.lossNN(x_train, y_train)\n",
        "        \n",
        "        loss.backward()\n",
        "                \n",
        "        self.iter += 1\n",
        "        \n",
        "        if self.iter % 100 == 0:\n",
        "        \n",
        "            print(loss)\n",
        "\n",
        "        return loss   "
      ],
      "metadata": {
        "id": "7EZjlbVMi8w0"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Generate data"
      ],
      "metadata": {
        "id": "pOe9yRvxjakj"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# get the analytical solution over the full domain\n",
        "x = torch.linspace(min,max,total_points).view(-1,1) #prepare to NN\n",
        "y = f_BC(x)\n",
        "print(x.shape, y.shape)\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Z1V03CRRjRML",
        "outputId": "debfb85f-3d11-47e0-faed-742169b53136"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "torch.Size([500, 1]) torch.Size([500, 1])\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "fig, ax1 = plt.subplots()\n",
        "ax1.plot(x.detach().numpy(),y.detach().numpy(),color='blue',label='Real_Train')\n",
        "#ax1.plot(x_train.detach().numpy(),yh.detach().numpy(),color='red',label='Pred_Train')\n",
        "ax1.set_xlabel('x',color='black')\n",
        "ax1.set_ylabel('f(x)',color='black')\n",
        "ax1.tick_params(axis='y', color='black')\n",
        "ax1.legend(loc = 'upper left')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 297
        },
        "id": "OE4RUZ27p0fT",
        "outputId": "4334362a-3564-4b8f-c44d-4f1b1f036a55"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7f752e39bc10>"
            ]
          },
          "metadata": {},
          "execution_count": 14
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZzVc///8cerqbRpH6G0IERfTRrhcv0IRS4SF1FCosVSZC1CpChFRCJFcUWI1GVLsl2WaIb2pKnIpDRFWVqnXr8/3p84TdM0c+ac8z7L6367nduc8zmfc87zWOY17897E1XFGGOMKakyvgMYY4xJTFZAjDHGhMUKiDHGmLBYATHGGBMWKyDGGGPCUtZ3gFiqXbu2NmzY0HcMY4xJKNnZ2etUNb3g8ZQqIA0bNiQrK8t3DGOMSSgi8kNhx+0SljHGmLBYATHGGBMWKyDGGGPCklJ9IIXZvn07ubm5bNmyxXeUpFGhQgXq1atHuXLlfEcxxkRRyheQ3Nxc9t9/fxo2bIiI+I6T8FSV9evXk5ubS6NGjXzHMcZEkddLWCLyrIisFZEFe3leRGSkiOSIyDwROS7kuS4isjS4dQk3w5YtW6hVq5YVjwgREWrVqmUtOmNSgO8+kPFA2yKePxtoHNx6AKMBRKQmMAA4AWgJDBCRGuGGsOIRWfbP05jU4PUSlqp+IiINizilPfC8ujXnZ4lIdRE5CGgFzFDVXwBEZAauEL0U3cSmMDt2wObNsHWru5+fDxs3wiOPQI0acOihcNhhULcuWG0xJnnEex9IXeDHkMe5wbG9Hd+DiPTAtV6oX79+dFKmmJ074bff/r4VdrVqwwa45Zbdjx18MPy//wdt2sD550OtWrHJa4yJDt+XsKJOVceoaqaqZqan7zETPy6kpaWRkZFB06ZNadeuHRs2bAjrfcaPH0+vXr0KfW7w4MFkZGSQkZHx1+dlZGQwcuTIYr13t27dyMpaxA8/wNy5kJMD69bBfvu5lsXhh8Mxx0CzZnDccVC/vmuF5OTAe+/ByJFw6qnw2WfQrRsceCD861/w1luuIBljEk+8t0BWAYeEPK4XHFuFu4wVevyjmKWKsIoVKzJnzhwAunTpwqhRo+jfv39EP6N///5/vWeVKlX++rxdVBVVpUyZPf+m+P13uP32sfz2G6xfD9WrQ+3aUKUKFHI64C5VVa3qbocd5lodvXuDKnzzDbzyCvznP3DuuXDEEdC3L1xxBZSN9/8ijTF/iff/XacBvURkEq7DfKOqrhaR6cADIR3nZwJ3lPbD+vSBAr9XSy0jAx59tPjnn3TSScybNw+AZcuWcf3115OXl0elSpV45plnOOqoo/jvf//LoEGD2LZtG7Vq1WLixInUqVOnxNm+//57zjrrLE444QSys7N5++23GTJkCLNnz2bz5s2cf/5FXHnlfWzcCNdc04pBg4bTunUm1atX4cYbb+TNN9+kYsWKTJ06tdifL+JaKMcdB/ffD5Mnw8MPw9VXu59DhriiYn0lxsQ/38N4XwK+AI4UkVwRuVpErhGRa4JT3gaWAznAM8B1AEHn+f3A7OA2cFeHeiLbsWMHM2fO5LzzzgOgR48ePP7442RnZzN8+HCuu+46AP75z38ya9YsvvnmGzp27MhDDz0U9mcuXbqU6667joULF9KgQQMGDx7MV19l8fbb83j33Y/5+ut51KvnWhu1a7sWwp9//smJJ57I3LlzOeWUU3jmmWfC+uxy5aBTJ5g9G157zXW+n3ceXHABrFoV9lcyxsSI71FYnfbxvALX7+W5Z4FnI5mnJC2FSNq8eTMZGRmsWrWKJk2a0KZNG/744w8+//xzOnTo8Nd5W7duBdzkx0suuYTVq1ezbdu2Uk3Ya9CgASeeeOJfj1944RWefnoM27fn88svq9m+fREHHnjsbq8pX7485557LgAtWrRgxowZYX8+uNbGv/8N7drBY4/BPfdAkybu/pVXWmvEmHiV9J3oiWBXH8gPP/yAqjJq1Ch27txJ9erVmTNnzl+3xYsXA9C7d2969erF/Pnzefrpp0s1aa9y5cqA65v46qsVDB8+nDFjZpKVNY927c4hP3/P9y5Xrtxfcz3S0tLIz88P+/N3f1+49VaYPx9atICrroLOnd1IL2NM/LECEkcqVarEyJEjefjhh6lUqRKNGjXi1VdfBVwn99y5cwHYuHEjdeu6UcsTJkwo9efm58OyZZCT8xtVqlSmZctqbNnyM++8806p3zschx0G778Pgwa5zvYWLSConcaYOGIFJM40b96cY489lpdeeomJEycybtw4mjVrxjHHHMPUqVMBuPfee+nQoQMtWrSgdu3apfo8Vfj2Wzfk9owzmnHCCc1p2vQoLr30Uk4++eRIfKWwpKVB//7w0UeuBXLiifDuu97iGGMKIa6bITVkZmZqwR0JFy9eTJMmTTwl8uvPP908jZ073V/9VatG7r0j+c915UrXuT5/PoweDT16RORtjTHFJCLZqppZ8Li1QFLUxo2wZInroD7qqMgWj0irXx8+/RTOPht69oShQ30nMsZA/M8DMWEYPHjwX30nu3To0OGviYQbNrg+jwoVoHFjKF/eR8qSqVIFpkxxo7L69YNff4UHH7QRWsb4ZAUE10GdTCvIhs46L2hX8ahY0c0Aj8bM72hdFi1XDl54AapVc62QLVtgxAgrIsb4kvIFpEKFCqxfvz4l9gTZVTwqVXItj2gVj/Xr11OhQoXIvzlu6ZRRo9waXI8+CpUrw+DBUfkoY8w+pHwBqVevHrm5ueTl5fmOElVbtsDate6v+AoVYOnS6H3Wri1to0XELRX/55/wwAOuiNx5Z9Q+zhizFylfQMqVK5f0W6/OnQutWrlVc//3v+RYRl3EjcjatMkN961eHYKVXowxMZLyBSTZLV8Obdu6UVbTpydH8dglLQ3Gj3fzRHr3dqO1ghVWjDExYMN4k9iGDW7PjW3bXPE45JB9vybRlC0LL70EzZvDJZfA11/7TmRM6rACkqTy890v1OXL3fDXo4/2nSh6KleGN9+E9HTXAlm50nciY1KDFZAkdeutbifA0aPhlFN8p4m+Aw90uxtu2uS2y9282XciY5KfFZAkNHasWwq9Tx+3UVOqOOYYmDjRbQrWo4db58sYEz1WQJJMVhZcfz2ceSYMG+Y7Teydcw7cd5/bLvfxx32nMSa5+d6RsK2ILBGRHBHpV8jzI0RkTnD7TkQ2hDy3I+S5abFNHp9+/RU6dIA6deDFF1N3f/H+/aF9e7j5ZvjkE99pjEle3n7FiEgaMApoA+QCs0Vkmqou2nWOqt4Ucn5voHnIW2xW1YxY5Y13qtC1K+TmJs9cj3CVKQPPPw/HH++2zJ07123Ha4yJLJ8tkJZAjqouV9VtwCSgfRHndwJeikmyBPTIIzB1qrtsFbJDbcqqWhVefhnWrXMLMFp/iDGR57OA1AV+DHmcGxzbg4g0ABoBH4QcriAiWSIyS0TO39uHiEiP4LysZF2u5IsvoG9ft6/4jTf6ThM/MjLg4Yfd6KzHHvOdxpjkkyid6B2Byaq6I+RYg2CDk0uBR0XksMJeqKpjVDVTVTPT09NjkTWmfv8dLrvMzcJ+9llbmbag6693/SG33w7Z2b7TGJNcfBaQVUDo3Oh6wbHCdKTA5StVXRX8XA58xO79IymjTx/4/vu/lzk3uxOBcePcwIKOHV3BNcZEhs8CMhtoLCKNRKQ8rkjsMZpKRI4CagBfhByrISL7BfdrAycDiwq+Ntm98YZrdfTrBx63L497tWq5+SHLlsFtt/lOY0zy8FZAVDUf6AVMBxYDr6jqQhEZKCLnhZzaEZiku+9S1ATIEpG5wIfAkNDRW6lgzRro3h2OOw4GDPCdJv6dcgrccgs8/bRbF8wYU3oSrd3j4lFmZqZmZWX5jlFqqm7Npw8+cIsHNmniO1Fi2LLFFdzffoMFC9wS8MaYfROR7KDPeTeJ0oluQowfD2+/DQ89ZMWjJCpUcPND1qyx0WrGRIIVkASzerWbYX3KKW6EkSmZzEy3e+Hzz7t5M8aY8FkBSSCqbte9LVvcgoll7N9eWO66y80R6dnTLf9ijAmP/QpKIJMnu5FXAwdC48a+0ySu8uXd6LV169z8EGNMeKyAJIj166FXL2jRAm66ad/nm6I1b+4uBY4dCx9/7DuNMYnJCkiCuOkm+OUX95dzqq6yG2kDBkCjRu5S1pYtvtMYk3isgCSAjz5yM8379oVjj/WdJnlUrux2bFyyBB54wHcaYxKPFZA4t22b6zhv1Mjtc2Ei66yzoHNnGDIEFi70ncaYxGIFJM6NGAGLF8PIkVCxou80yemRR2D//V2hTqF5tcaUmhWQOLZypRtxdf75bua5iY4DDnCXsD75BF6yHWeMKTYrIHGsTx/389FH/eZIBd26uRFut95qK/YaU1xWQOLUW2/BlClwzz3QoIHvNMkvLQ1GjXIz/e+/33caYxKDFZA4tGUL9O7t1rmyOR+xc8IJcNVVf/c7GWOKZgUkDo0YAStWwBNPuFnTJnaGDIEqVeCGG6xD3Zh9sQISZ1avhsGDXcf56af7TpN60tPdJaz334fXXvOdxpj4ZgUkzvTv7+Z+DB/uO0nquuYaaNbMbUC1ebPvNMbEL68FRETaisgSEckRkX6FPH+liOSJyJzg1i3kuS4isjS4dYlt8ujIznZ7ffTpA4cd5jtN6ipb1o18W7nSRsAZUxRvOxKKSBrwHdAGyMXtkd4pdGtaEbkSyFTVXgVeWxPIAjIBBbKBFqpa5OLc8bwjoarb4+O772DpUqha1Xcic/75MHMm5ORAnTq+0xjjTzzuSNgSyFHV5aq6DZgEtC/ma88CZqjqL0HRmAG0jVLOmHj1Vfj0Uxg0yIpHvHjoITci7p57fCcxJj75LCB1gR9DHucGxwq6UETmichkETmkhK9FRHqISJaIZOXl5UUid8Rt3uz2pWjWzA0jNfHhiCPcro9jx7o91I0xu4v3TvT/Ag1V9VhcK2NCSd9AVceoaqaqZqanp0c8YCQ88gj88IMbvpuW5juNCXXPPVCtmpuhbozZnc8Csgo4JORxveDYX1R1vapuDR6OBVoU97WJ4uef4cEH3fX2007zncYUVLOmKyLTp8O77/pOY0x88VlAZgONRaSRiJQHOgLTQk8QkYNCHp4H7JofPB04U0RqiEgN4MzgWMIZONBdZx861HcSszfXXQeHH+6G9ebn+05jTPzwVkBUNR/ohfvFvxh4RVUXishAETkvOO0GEVkoInOBG4Arg9f+AtyPK0KzgYHBsYSydCmMGQM9erjr7SY+lS8Pw4bBokXwzDO+0xgTP7wN4/Uh3obxXnKJWzQxJwcOPNB3GlMUVWjVCr79FpYtc8udGJMq4nEYb0qbPRteecVdFrHiEf9E3GXGtWttcqExu1gB8UDVDdtNT3cFxCSGE090gx2GDYN163ynMcY/KyAeTJ8OH30Ed99tkwYTzeDB8McfbuScManOCkiM7dwJffvCoYdCz56+05iSOvpo6NLFbT61cqXvNMb4ZQUkxiZOhHnz3F+yttdHYrr3XncZ8t57fScxxi8rIDG0ZQvcdZfbe/vii32nMeGqX98tcTJhghvaa0yqsgISQ0895S57DB0KZeyffEK7806oXNnt32JMqrJfYzGyq+P1jDPczSS22rXhttvgjTdg1izfaYzxwwpIjDzxhJtDcP/9vpOYSLnpJjjgAOjXz/ZPN6nJCkgMbNzo9pY45xw46STfaUykVKnihmJ//DG8957vNMbEnhWQGHj0Ufj1V7dwokku3bu7TvW777ZWiEk9VkCibP16t9/HhRfCccf5TmMibb/9XPGYPduta2ZMKrECEmXDh8Pvv8N99/lOYqKlSxc3MfSee6wVYlKLFZAo+vlnGDkSOnWCY47xncZES7lyrnh88w1MmeI7jTGxYwUkioYMga1bYcAA30lMtHXu7PZ0GTDALVdjTCrwWkBEpK2ILBGRHBHpV8jzN4vIIhGZJyIzRaRByHM7RGROcJtW8LW+5ebC6NHu8oZtFpX8ypZ1S5ssWACvvuo7jTGx4W1DKRFJA74D2gC5uJ0FO6nqopBzTgO+VNVNInIt0EpVLwme+0NVS7StTyw3lLr2Whg3Dr77Dho2jMlHGs927oRjj4UdO1whSUvznciYyIjHDaVaAjmqulxVtwGTgPahJ6jqh6q6KXg4C6gX44xhWbECxo6Fbt2seKSSMmXcYIlvv4UXX/Sdxpjo81lA6gI/hjzODY7tzdXAOyGPK4hIlojMEpHz9/YiEekRnJeVl5dXusTFdP/97q9PWycp9VxwAWRkuEKyfbvvNMZEV0J0oovIZUAmMCzkcIOgSXUp8KiIHFbYa1V1jKpmqmpmenp61LMuXepWab32WqhbVDk0SWlXK2TZMnjhBd9pjIkunwVkFXBIyON6wbHdiEhroD9wnqpu3XVcVVcFP5cDHwHNoxm2uAYPdpPL+u0xJMCkinbt4Pjj3coD27b5TmNM9PgsILOBxiLSSETKAx2B3UZTiUhz4Glc8VgbcryGiOwX3K8NnAx435lh2TL4z3/gmmugTh3faYwvIq54/PADPPec7zTGRI+3AqKq+UAvYDqwGHhFVReKyEAROS84bRhQBXi1wHDdJkCWiMwFPgSGhI7e8mXwYDep7LbbfCcxvp11FpxwglvC31ohJll5G8brQzSH8a5YAY0bu53qHnssKh9hEsw778C//gXPPONG5BmTqOJxGG9SeeABN5msb1/fSUy8aNvW9YU88ICNyDLJyQpIBPzwA4wf7/7KPPhg32lMvBBxa2StWOH6xoxJNlZAIuDBB93wTRt5ZQo65xxo0cL1j+Xn+05jTGRZASmllSvh2Wfh6quhXkLMkzextKsVsmwZTJzoO40xkWUFpJSGDnU/rfVh9qZdOzc73VohJtlYASmF3Fy35lXXrm5bU2MKs6sVsnQpTJrkO40xkWMFpBSGDnUrsN5xh+8kJt61b+9W6h00yK3Wa0wysAISpp9+cuP7u3SxFXfNvpUp41ohS5bAK6/4TmNMZFgBCdNDD7nr2Xfe6TuJSRQXXABNm7rVmq0VYpKBFZAwrF4NTz8Nl18Ohx7qO41JFGXKwN13w+LFMHmy7zTGlJ4VkDAMH+5mFtt+H6akLroIjj7atUJs73ST6KyAlNDPP7u9zjt3hsMP953GJJoyZeCuu2DhQnj9dd9pjCkdKyAlNHw4bN1qrQ8TvosvhqOOcku+WyvEJDIrICWQlwdPPgmdOsERR/hOYxJVWpprhcyfD2+84TuNMeErVgERkQNE5AIRuV5ErhKRliKScsXn4Ydh82b3P78xpXHJJW75/4EDIYV2VDBJpsgiICKnich04C3gbOAg4GjgLmC+iNwnIlWjH9O/devgiSfc//hHHeU7jUl0Zcu6P0TmzoVp0/Z9vjHxaF+tiH8B3VX1eFXtoap3qeqtqnoe0Az4BmgT7oeLSFsRWSIiOSKyx2pSIrKfiLwcPP+liDQMee6O4PgSETkr3AzFNWIEbNrkhmEaEwmXXgqHHWatEJO4iiwgqnqbqq7cy3P5qvqGqr4WzgeLSBowCteyORroJCJHFzjtauBXVT0cGAEMDV57NG4P9WOAtsCTwftFxS+/wOOPQ4cObgimMZFQtqwbjPH11/DWW77TGFNyxe0DeUFEqoU8bigiM0v52S2BHFVdrqrbgElA+wLntAcmBPcnA2eIiATHJ6nqVlVdAeQE7xcVI0bA779b34eJvMsug0aNrBVioufbb92+NMuXR/69i9sR/inwpYj8S0S6A+8Bj5bys+sCP4Y8zg2OFXqOquYDG4FaxXwtACLSQ0SyRCQrLy8vrKDr1rmhl//3f2G93Ji9KlfOLYczeza8+67vNCYZDR4MH30E++8f+fcuVgFR1aeBbsBUYCBwiqr+N/JxIk9Vx6hqpqpmpqenh/Ueo0fDiy9GOJgxgSuugAYN4L77rBViImvpUve767rrIMxff0Uq7iWsy4FngSuA8cDbItKslJ+9Cjgk5HG94Fih54hIWaAasL6Yr42otKj1sJhUV7682xLgyy9hxgzfaUwyefBB99/XLbdE5/2LewnrQuCfqvqSqt4BXIMrJKUxG2gsIo1EpDyuU7zggMZpQJfg/kXAB6qqwfGOwSitRkBj4KtS5jHGm65d4ZBDrBViImfFCnj+eejZEw48MDqfUdxLWOer6tqQx18BJ5Tmg4M+jV7AdGAx8IqqLhSRgSJyXnDaOKCWiOQANwP9gtcuBF4BFgHvAterqi2QbRLWrlbI55/DBx/4TmOSwYMPupF+t98evc8QLeLPHRG5C3hSVX/Zy/OnA5VU9c0o5YuozMxMzcrK8h3DmEJt3ermhRx6KHz8sdsK15hwrFzpFnvt3h1GjSr9+4lItqpmFjxedh+vmw/8V0S2AF8DeUAF3CWjDOB94IHSxzPG7Lcf9OsHvXu7AtKqle9EJlENHep+9u0b3c/Z1yWsi1T1ZNxlpoVAGvAb8B+gparepKrhjY01xuyhWzc46CDXF2JMOFatgnHj4MoroX796H7WvlogLUTkYKAzcFqB5yoCm6OSypgUVaGC+6uxTx/45BM45RTfiUyiGTrUbZkci+2299UCeQqYCRwFZIXcsoOfxpgI69ED6tRxs9ONKYmffoIxY6BLF2jYMPqft6+1sEaqahPgWVU9NOTWSFVtN3BjoqBiRTdyZuZM+Owz32lMIhk6FPLzY9P6gOIP47022kGMMX/r2dPNHLZWiCmu1atd6+OKK9xIvlhIuU2hjEkElSvDbbfBe+/BrFm+05hE8NBDsH17bLfbtgJiTJy69lqoXdtaIWbf1qyBp56Cyy93c4lixQqIMXGqShW3htE778BXtlCPKcKwYbFvfYAVEGPi2vXXQ82a1goxe/fzz27F8M6d3ezzWLICYkwc239/uPlmt2NhdrbvNCYeDRvmlsHxseGdFRBj4lzv3lC9urVCzJ7WroUnn4RLL4XGjWP/+VZAjIlzVavCTTfBtGnwzTe+05h4Mny4v9YHWAExJiHccANUqwb33+87iYkXeXlupd2OHeHII/1ksAJiTAKoXh1uvBGmTIF583ynMfHg4Ydh82Z/rQ+wAmJMwujTx13OslaIWbcOnnjCtT6aNPGXw0sBEZGaIjJDRJYGP2sUck6GiHwhIgtFZJ6IXBLy3HgRWSEic4JbRmy/gTGxV6OGu5Q1eTIsWOA7jfHpkUdg0ya4+26/OXy1QPoBM1W1MW61336FnLMJuEJVjwHaAo+KSPWQ529T1YzgNif6kY3xr08fN8HQWiGpa906ePxxuPhiv60P8FdA2gMTgvsTgPMLnqCq36nq0uD+T8BaID1mCY2JQ7VquVbIq6/C/Pm+0xgfHnrItT4GDPCdxF8BqaOqq4P7a4A6RZ0sIi2B8sCykMODg0tbI0RkvyJe20NEskQkKy/PNk80ie+WW9wEw3j4BWJia80a1/fRubP/1gdEsYCIyPsisqCQW/vQ81RVAS3ifQ4CXgC6qurO4PAduE2ujgdqAnvd+VdVx6hqpqpmpqdbA8Ykvpo1XRGZMsVmp6eaBx+Ebdvgnnt8J3GiVkBUtbWqNi3kNhX4OSgMuwrE2sLeQ0SqAm8B/VV1Vsh7r1ZnK/Ac0DJa38OYeNSnjyskvjtRTez8+KNbcbdr19ivebU3vi5hTQO6BPe7AFMLniAi5YEpwPOqOrnAc7uKj+D6T2xMikkpVau6XQvfeQc+/9x3GhMLgweDqt95HwX5KiBDgDYishRoHTxGRDJFZGxwzsXAKcCVhQzXnSgi84H5QG1gUGzjG+Nfr15wwAHWCkkFy5fDuHHQowc0aOA7zd/EdUGkhszMTM3KyvIdw5iIeewxdzlr5kw4/XTfaUy0dO0KkybBsmVw8MGx/3wRyVbVzILHbSa6MQmsZ0+oW9e1QlLob8GUsmQJPP88XHedn+JRFCsgxiSwChXcNfHPP4fp032nMdFw331QsSL03etYU3+sgBiT4K66Cho2dIXEWiHJZcECd+nqhhtcf1e8sQJiTIIrX95NKszOhql7jGc0iWzAADdp9NZbfScpnBUQY5LAZZfBEUe4vpAdO3ynMZEweza8/rrb0rhmTd9pCmcFxJgkULasu1a+YAG8+KLvNKa0VKFfP0hPdwUkXlkBMSZJXHwxHHec6wvZssV3GlMaM2bABx+4f5f77+87zd5ZATEmSZQpA0OHwsqVMHq07zQmXDt3uhFXDRu6YdrxzAqIMUmkdWto0wYGDYKNG32nMeF4+WWYM8ft+bLfXtcZjw9WQIxJMkOGwC+/uH0jTGLZts1dtmrWDC691HeafbMCYkySOe446NQJRoyAn37yncaUxJgxbt2rBx90lyTjXQJENMaU1KBBkJ/vRmaZxPDHH+6y1amnQtu2vtMUjxUQY5LQoYfCNde4FVy//dZ3GlMcjzwCa9e6S5AivtMUjxUQY5LUXXe5NZT69/edxOzLmjUwbBhccAGceKLvNMVnBcSYJHXAAXDbbW428xdf+E5jinLPPW7uztChvpOUjJcCIiI1RWSGiCwNftbYy3k7QjaTmhZyvJGIfCkiOSLycrB7oTGmgJtvhjp13B7qttBifJo3z11q7NULGjf2naZkfLVA+gEzVbUxMDN4XJjNqpoR3M4LOT4UGKGqhwO/AldHN64xialKFXjgAdcCefll32lMQaquuFerlpg7S/oqIO2BCcH9Cbh9zYsl2Af9dGDXPukler0xqaZLF2je3O2hvnmz7zQm1Ntvw/vvu1V343XBxKL4KiB1VHV1cH8NUGcv51UQkSwRmSUiu4pELWCDquYHj3OBunv7IBHpEbxHVl5eXkTCG5NI0tLcnJAff4SHH/adxuyyfbtbpv2II9xug4mobLTeWETeBw4s5KndxoSoqorI3q7ONlDVVSJyKPCBiMwHSrRAg6qOAcaA2xO9JK81JlmceipceKGboHbVVfG3NWoqGjPGDbGeOhXKlfOdJjxRa4GoamtVbVrIbSrws4gcBBD8XLuX91gV/FwOfAQ0B9YD1UVkV/GrB6yK1vcwJlk89JCbXHjnnb6TmA0b3GWr02nR1scAAA4ZSURBVE6Ddu18pwmfr0tY04Auwf0uwB77qIlIDRHZL7hfGzgZWKSqCnwIXFTU640xuzv0ULjpJpgwAbKyfKdJbYMGufXKHnkkcSYNFsZXARkCtBGRpUDr4DEikikiY4NzmgBZIjIXVzCGqOqi4Lm+wM0ikoPrExkX0/TGJKg773TzQ/r0sWG9vixeDI89Bl27QkaG7zSlI5pC/xVlZmZqlv3pZVLc2LHQvbvbubBTJ99pUouqW24/OxuWLHHFPBGISLaqZhY8bjPRjUkxXbtCZqabf/Dbb77TpJbJk2HmTHcJK1GKR1GsgBiTYtLS4Mkn3fpL997rO03q+OMPtzJARoZb6DIZWAExJgUdfzz06AEjR7qlNEz0DRoEubkwapQr4snACogxKeqBB6BGDTeJbedO32mS25IlbsTVlVfCP/7hO03kWAExJkXVrOnmhnz2GTz/vO80yUsVeveGSpUSb7XdfbECYkwK69LF/UV8++1uXoKJvJdeghkz3G6DydBxHsoKiDEprEwZ16G+fj3ccYfvNMln/Xo356Zly8Rd76ooVkCMSXHNmrkZ6mPGwEcf+U6TXG65BX79FZ55Jnk6zkNZATHGMHCgW+qke3db8j1S3n/fLRtz++1w7LG+00SHFRBjDJUquRnqOTk2NyQSNm2Cnj3dDoOJuFFUcVkBMcYAbmXYbt1g+HBbbLG07r0Xli93lwUrVPCdJnqsgBhj/jJsmNtD/eqr3YZHpuQ+/9xt3NW9O7Rq5TtNdFkBMcb8pXp1GD3azU4fNMh3msTz559uaHT9+qmx+6MVEGPMbtq3hyuugMGD4csvfadJLP36uX6k556D/ff3nSb6rIAYY/YwciTUrQuXX+7+qjb7NnMmPPGEm/eR7JeudrECYozZQ7VqbnmTnBy47TbfaeLfxo1umfwjj3RrjKUKLwVERGqKyAwRWRr8rFHIOaeJyJyQ2xYROT94bryIrAh5LsH39TIm/px6qpsIN3o0vPOO7zTxS9Utz/7TT67oVqzoO1Hs+GqB9ANmqmpjYGbweDeq+qGqZqhqBnA6sAl4L+SU23Y9r6pzYpLamBRz//3QtKn763rNGt9p4tOzz8KkSW4yZsuWvtPElq8C0h6YENyfAJy/j/MvAt5R1U1RTWWM2U2FCm4xwN9+g86dYccO34niy8KFbqXd1q1dB3qq8VVA6qjq6uD+GqDOPs7vCLxU4NhgEZknIiNEZL+9vVBEeohIlohk5eXllSKyMampaVO3CdIHH7gWiXE2bYJLLnGjrV54wS1MmWqi9pVF5H0RWVDIrX3oeaqqgBbxPgcB/wdMDzl8B3AUcDxQE+i7t9er6hhVzVTVzPT09NJ8JWNS1pVXuqG9Awe6NZ6MW4By4UJXPA480HcaP8pG641VtfXenhORn0XkIFVdHRSItUW81cXAFFX9a15sSOtlq4g8B9wakdDGmEKJuGXfs7Lcpaw5c+Cgg3yn8mf8eLdMSb9+cOaZvtP446vRNQ3oEtzvAkwt4txOFLh8FRQdRERw/ScLopDRGBOicmV49VX44w+46CLYutV3Ij+++sqNujrjDLuk56uADAHaiMhSoHXwGBHJFJGxu04SkYbAIcDHBV4/UUTmA/OB2oAtumBMDBx9tJtl/fnnboMk3evF5+S0Zg38+9+u9fXyy1A2atdwEoOXr6+q64EzCjmeBXQLefw9ULeQ806PZj5jzN5dfDHMn+/Wyjr2WLjxRt+JYmPrVujQwW39+8UXUKuW70T+pXj9NMaE4777YMECuPlmaNIk+fsBVOGqq+DTT92cj2bNfCeKDyk48MwYU1plyrjRR02buv6Qb77xnSi67roLXnzRLVNyySW+08QPKyDGmLBUqQJvvw01asDZZ7sNlJLRmDGucHTvnpqTBYtiBcQYE7a6deHdd93mU2edBWuLGpCfgF5/3Q0WaNvWDWMW8Z0ovlgBMcaUSpMm8OabsGqVKyLr1/tOFBlvvQUdO7r1rV591UZcFcYKiDGm1E46yf21vngxtGnjRiolshkz4MILXWf5O++4y3VmT1ZAjDER0bYtvPEGLFrkFhdM1CIyfbrblfHII939atV8J4pfVkCMMRETWkRatXKXtRLJyy9Du3ZwxBGuFVKzpu9E8c0KiDEmotq2dX0iK1a4S1uLFvlOVDyjR0OnTnDiifDxx3DAAb4TxT8rIMaYiGvdGj75xI3OOvlk9ws5Xu3YAbff7kZbnXuuXbYqCSsgxpioaN7cLflx4IGuoDz+ePytnbVhgysaw4a5AvLaa6m1JW1pWQExxkRNw4auiJx9NtxwA1x2Gfz+u+9UzpdfQosWbn+Tp592m2aVK+c7VWKxAmKMiarq1V3H+uDBf68j9ckn/vLk58OQIfDPf7rLVx9/DD16+MuTyKyAGGOirkwZuPNOVzjKlHEjtG680V1CiqWsLDcx8I473LLsc+bAP/4R2wzJxAqIMSZmTj7Z/dK+7jrXJ3LEETB2rGsVRFNurmtlnHCC29Nj8mTXGqpePbqfm+ysgBhjYqpKFXjiCcjOdgWke3c3aW/sWNi2LbKf9f33bsn5ww9329D26uVmy194oa1rFQleCoiIdBCRhSKyU0QyizivrYgsEZEcEekXcryRiHwZHH9ZRMrHJrkxJlKaN4f//c/1j9So4QpJ3bpwyy1uw6pwR2xt3gxTp7pLVIcdBiNHujWtvvsOHnvMhuhGkqiHcXUi0gTYCTwN3BrsRFjwnDTgO6ANkAvMBjqp6iIReQV4XVUnichTwFxVHb2vz83MzNSsrD0+yhjjmaqb+T1mjPvln58PDRq40Vv/+IfreD/ySNhvvz1fl5cHS5fCrFnw2Wfw3nvw559Qu7YrStdeC4cc4ud7JQsRyVbVPf7Y97Wl7WIAKboN2RLIUdXlwbmTgPYishg4Hbg0OG8CcC+wzwJijIlPIm5XwzPPdEvCT5niFjF84QV46qm/z6ta9e+tZLdvdyv/bt789/ONGkHnzm6Tq1atbFhutMXzAsV1gR9DHucCJwC1gA2qmh9yfI9903cRkR5AD4D69etHJ6kxJmIOOAB69nS3/Hx36WnOHFi2zLU21q93I7nKlXOXvho0cIXj+OPdpEUTO1ErICLyPlDYv87+qjo1Wp9bkKqOAcaAu4QVq881xpRe2bJw9NHuZuJP1AqIqrYu5VusAkKvXNYLjq0HqotI2aAVsuu4McaYGIrnYbyzgcbBiKvyQEdgmrpe/w+Bi4LzugAxa9EYY4xxfA3jvUBEcoGTgLdEZHpw/GAReRsgaF30AqYDi4FXVHVh8BZ9gZtFJAfXJzIu1t/BGGNSnZdhvL7YMF5jjCm5vQ3jjedLWMYYY+KYFRBjjDFhsQJijDEmLFZAjDHGhCWlOtFFJA/4IcyX1wbWRTBOrCV6fkj875Do+SHxv0Oi5wc/36GBqqYXPJhSBaQ0RCSrsFEIiSLR80Pif4dEzw+J/x0SPT/E13ewS1jGGGPCYgXEGGNMWKyAFN8Y3wFKKdHzQ+J/h0TPD4n/HRI9P8TRd7A+EGOMMWGxFogxxpiwWAExxhgTFisgxSAibUVkiYjkiEg/33lKQkSeFZG1IrLAd5ZwiMghIvKhiCwSkYUicqPvTCUlIhVE5CsRmRt8h/t8ZwqHiKSJyDci8qbvLOEQke9FZL6IzBGRhFtVVUSqi8hkEflWRBaLyEneM1kfSNFEJA34DmiD2z53NtBJVRd5DVZMInIK8AfwvKo29Z2npETkIOAgVf1aRPYHsoHzE+WfP4CICFBZVf8QkXLAp8CNqjrLc7QSEZGbgUygqqqe6ztPSYnI90CmqibkREIRmQD8T1XHBnskVVLVDT4zWQtk31oCOaq6XFW3AZOA9p4zFZuqfgL84jtHuFR1tap+Hdz/Hbc3TF2/qUpGnT+Ch+WCW0L95SYi9YBzgLG+s6QiEakGnEKw95GqbvNdPMAKSHHUBX4MeZxLgv0CSxYi0hBoDnzpN0nJBZd/5gBrgRmqmmjf4VHgdmCn7yCloMB7IpItIj18hymhRkAe8FxwGXGsiFT2HcoKiEkIIlIFeA3oo6q/+c5TUqq6Q1UzgHpASxFJmMuJInIusFZVs31nKaV/qupxwNnA9cHl3URRFjgOGK2qzYE/Ae/9sVZA9m0VcEjI43rBMRMjQb/Ba8BEVX3dd57SCC47fAi09Z2lBE4Gzgv6ECYBp4vIf/xGKjlVXRX8XAtMwV2eThS5QG5Iy3UyrqB4ZQVk32YDjUWkUdBx1RGY5jlTygg6oMcBi1X1Ed95wiEi6SJSPbhfETcg41u/qYpPVe9Q1Xqq2hD33/8HqnqZ51glIiKVg0EYBJd+zgQSZmSiqq4BfhSRI4NDZwDeB5KU9R0g3qlqvoj0AqYDacCzqrrQc6xiE5GXgFZAbRHJBQao6ji/qUrkZOByYH7QhwBwp6q+7TFTSR0ETAhG9JUBXlHVhBwKm8DqAFPc3yOUBV5U1Xf9Riqx3sDE4A/Z5UBXz3lsGK8xxpjw2CUsY4wxYbECYowxJixWQIwxxoTFCogxxpiwWAExxhgTFisgxhhjwmIFxBhjTFisgBjjkYgcLyLzgj1DKgf7hSTMOlkmtdlEQmM8E5FBQAWgIm69owc9RzKmWKyAGONZsDTFbGAL8A9V3eE5kjHFYpewjPGvFlAF2B/XEjEmIVgLxBjPRGQabpn0Rrjte3t5jmRMsdhqvMZ4JCJXANtV9cVgtd7PReR0Vf3AdzZj9sVaIMYYY8JifSDGGGPCYgXEGGNMWKyAGGOMCYsVEGOMMWGxAmKMMSYsVkCMMcaExQqIMcaYsPx/sZGQfSg2y9kAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#def get_training_data(x):\n",
        "#Nu: Number of training point, # Nf: Number of colloction points\n",
        "# Set Boundary conditions:\n",
        "BC_1=x[0,:]\n",
        "BC_2=x[-1,:]\n",
        "# Total Training points BC1+BC2\n",
        "all_train=torch.vstack([BC_1,BC_2])\n",
        "#Select Nu points\n",
        "idx = np.random.choice(all_train.shape[0], Nu, replace=False) \n",
        "x_BC=all_train[idx]\n",
        "#Select Nf points\n",
        "# Latin Hypercube sampling for collocation points \n",
        "x_PDE = BC_1 + (BC_2-BC_1)*lhs(1,Nf)\n",
        "x_PDE = torch.vstack((x_PDE,x_BC)) "
      ],
      "metadata": {
        "id": "_JMjSMOKdlWt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "x_PDE.shape"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SxRNKJ4Ejif8",
        "outputId": "dd5d17a4-00e3-4868-b2b5-c3a36efe1ccc"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "torch.Size([252, 1])"
            ]
          },
          "metadata": {},
          "execution_count": 16
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Train Neural Network"
      ],
      "metadata": {
        "id": "I82eH55ElVSq"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#Store tensors to GPU\n",
        "torch.manual_seed(123)\n",
        "x_PDE=x_PDE.float().to(device)\n",
        "x_BC=x_BC.to(device)\n",
        "#Create Model\n",
        "model = FCN(layers)\n",
        "print(model)\n",
        "model.to(device)\n",
        "params = list(model.parameters())\n",
        "optimizer = torch.optim.Adam(model.parameters(),lr=lr,amsgrad=False)\n",
        "start_time = time.time()"
      ],
      "metadata": {
        "id": "8KZvgKcalOVV",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "903bed15-1471-4dfa-b812-2dbabb307882"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "FCN(\n",
            "  (activation): Tanh()\n",
            "  (loss_function): MSELoss()\n",
            "  (linears): ModuleList(\n",
            "    (0): Linear(in_features=1, out_features=50, bias=True)\n",
            "    (1): Linear(in_features=50, out_features=50, bias=True)\n",
            "    (2): Linear(in_features=50, out_features=20, bias=True)\n",
            "    (3): Linear(in_features=20, out_features=50, bias=True)\n",
            "    (4): Linear(in_features=50, out_features=50, bias=True)\n",
            "    (5): Linear(in_features=50, out_features=1, bias=True)\n",
            "  )\n",
            ")\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "for i in range(steps):\n",
        "    yh = model(x_PDE)\n",
        "    loss = model.loss(x_PDE,x_BC)# use mean squared error\n",
        "    optimizer.zero_grad()\n",
        "    loss.backward()\n",
        "    optimizer.step()\n",
        "    if i%(steps/10)==0:\n",
        "      print(loss)"
      ],
      "metadata": {
        "id": "KM3GHnh9ljCh",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "89c09060-a4c8-4b76-abea-3f34e0f301c2"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "tensor(1.5134, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0001, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(8.4055e-05, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(0.0062, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(7.3839e-05, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(6.3905e-05, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(5.1710e-05, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(3.9399e-05, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(2.8976e-05, device='cuda:0', grad_fn=<AddBackward0>)\n",
            "tensor(2.1685e-05, device='cuda:0', grad_fn=<AddBackward0>)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Plots"
      ],
      "metadata": {
        "id": "kpRh89zM8Te8"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Function\n",
        "yh=model(x.to(device))\n",
        "y=f_BC(x)\n",
        "#Error\n",
        "print(model.lossBC(x.to(device)))"
      ],
      "metadata": {
        "id": "6kJ3gQ6WYhp9",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "6087df9b-625e-4435-a541-82c1ffb5341d"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "tensor(1.6701e-05, device='cuda:0', grad_fn=<MseLossBackward0>)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Derivative\n",
        "g=x.to(device)\n",
        "g=g.clone()\n",
        "g.requires_grad=True #Enable differentiation\n",
        "f=model(g)\n",
        "f_x=autograd.grad(f,g,torch.ones([g.shape[0],1]).to(device),retain_graph=True, create_graph=True)[0]"
      ],
      "metadata": {
        "id": "jEUUqEjr8Sim"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Detach from GPU\n",
        "y_plot=y.detach().numpy()\n",
        "yh_plot=yh.detach().cpu().numpy()\n",
        "f_x_plot=f_x.detach().cpu().numpy()"
      ],
      "metadata": {
        "id": "udWgFd6mqY7B"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# Plot\n",
        "fig, ax1 = plt.subplots()\n",
        "ax1.plot(x,y_plot,color='blue',label='Real')\n",
        "ax1.plot(x,yh_plot,color='red',label='Predicted')\n",
        "ax1.plot(x,f_x_plot,color='green',label='Derivative')\n",
        "ax1.set_xlabel('x',color='black')\n",
        "ax1.set_ylabel('f(x)',color='black')\n",
        "ax1.tick_params(axis='y', color='black')\n",
        "ax1.legend(loc = 'upper left')"
      ],
      "metadata": {
        "id": "M9YZDGcpodnQ",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 297
        },
        "outputId": "bd218024-30e2-4073-bf0f-837910d10185"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.legend.Legend at 0x7f752cb19510>"
            ]
          },
          "metadata": {},
          "execution_count": 22
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1RUx9/H8ffQqygIgmBBxYINFXuvMdVfokZNsbfYokZjjWLDXqPRGI0xJpaYptEYjVijsWDvghVQEbHQy+7O8wfoQ6yIwN2FeZ2zB/buLR/W43537p07I6SUKIqiKMrLMtM6gKIoimKaVAFRFEVRskQVEEVRFCVLVAFRFEVRskQVEEVRFCVLLLQOkJsKFy4sS5YsqXUMRVEUk3LkyJE7UkrXx5fnqwJSsmRJgoODtY6hKIpiUoQQ1562XJ3CUhRFUbJEFRBFURQlS1QBURRFUbIkX10DeZrU1FTCw8NJSkrSOkqeYmNjg5eXF5aWllpHURQlh+T7AhIeHo6joyMlS5ZECKF1nDxBSkl0dDTh4eF4e3trHUdRlByS709hJSUl4eLioopHNhJC4OLiolp1ipLHaVpAhBDfCiFuCyFOP+N1IYRYIIQIFUKcFEJUz/BaFyFESPqjyyvmeJXNladQ76mi5H1an8L6DlgIfP+M118HfNIftYHFQG0hhDMwHvAHJHBECLFRSnkvJ8NKKUk1pJKQkkCKIQWdQYdAYG5mjrW5NbaWtliZW+VkBEVRFKOhaQGRUu4RQpR8ziptgO9l2qQlB4QQBYUQHkAT4G8p5V0AIcTfQGtgTU7kjE+J527iXe4l3SNFn/LcdW0sbChoU5DCdoWxsbDJ1P7Nzc2pXLkyOp0Ob29vVq1aRcGCBV8653fffUdwcDALFy586W1fVWoqJCeDTgd6fdqyhATYtQtKlQJPTzA3z/VYipKvSCm5l3SPK/euEBEbQUJqAompiQgh+F/5/1HQ5uU/V55H6xbIi3gCYRmeh6cve9byJwghegO9AYoXL56lEBGxEcQmx1LAugBF7Itgb2WPtbk1FmYWSCQ6g45kXTLxqfHEJMdwK+4Wt+Ju4WTtRFHHothb2T93/7a2thw/fhyALl26sGjRIsaMGZOlrLlBSkhMlMTfTUb3IB6RnISVIQkLdFiixwqJxAyzO7cIez2QPXhz0aIiCVXqUKZZcVq0FDRtCqqDlqK8ukt3L7Hxwkb2Xt/L/rD9RMZHPnW9ul51810BeWVSyqXAUgB/f/8sTb9Y3Kk4FmYWWJg9+XYJBFbmVliZW+Fo7Yi7gzsp+hTuJNzhdvxtzt05h7OtM8UKFMPS/MWfmHXr1uXkyZMAXLp0if79+xMVFYWdnR3ffPMN5cuX548//mDy5MmkpKTg4uLCjz/+SJEiRbLyp72UlGRJ7I0YxP17OOhjcCWtNSYBvYU1WFiChSXCTIDBgFUMtHPbg3XUasx0BjgK4Uc92TTrLb5xfJfC7zejVz9Lqld//nEVRfmvmOQYvj32LatOruLozaMAlCpUilalW1G1SFVKFSpFMadiOFg5YGNhg5QSzwJP/Y79Soy9gEQAxTI890pfFkHaaayMy3e96sEGD4b0hsBjMncq6v9ZAUWReODjG88nYy/wIOkBxZ2K42Ln8syt9Ho9QUFB9OjRA4DevXuzZMkSfHx8OHjwIP369WPHjh00aNCAAwcOIIRg2bJlzJgxg9mzZ79kxsxLikkhKew29ol3cEGHHnN0dgXQFXLHwskBYWODhdmT/THMpB7byGuQkgKnT8O//+KxfSc9//qBvrFfE/ltEZYt78FU/150DSjJG2+AuvauKM8WGRfJ9H3TWXZ0GbEpsdQsWpNZLWfRzrcdJQqWyPU8xl5ANgIDhBBrSbuI/kBKeVMIsRUIFEIUSl+vFTBKq5DPIhA4WDlQ0a0iV+9f5cr9K8SmxFKsQDHMzf7/gkBiYiJ+fn5ERERQoUIFWrZsSVxcHPv376d9+/aP1ktOTgbS7l3p0KEDN2/eJCUlJcfutUiNTybp8k3sk6OxRpJoXRDcXbB0ccL8KQXjmaysoHp1qF4d8/79ISkJtm3DefEyRm2dhiF4Oqve+phONcYwcH4Z6tfPkT9HUUxWfEo8s/+dzcz9M0lMTaRDpQ4MqTME/6L+2gaTUmr2IO2i900glbTrGD2AvkDf9NcFsAi4BJwC/DNs2x0ITX90y8zxatSoIR939uzZJ5blBL1BL8MehMnDEYflmdtnZLIu+dFr9vb2Ukop4+PjZYMGDeT8+fPlgwcPpLu7+1P31bhxY7lhwwYppZQ7d+6UjRs3llJKuWLFCtm/f/9XzmpI1cnYc2FSfzhY6g8Hy9gz12RqXOJL7ydT7+21a1I3aLBMtbSRqZjLRXwi+3W4I6OishBcUfKgvy/9LUvOKykJQLZd11ZeuHMh1zMAwfIpn6ma3gcipewkpfSQUlpKKb2klMullEuklEvSX5dSyv5SytJSyspSyuAM234rpSyT/lih3V+ROWbCDK8CXpRxLkOSLonzd86TmJr4n3Xs7OxYsGABs2fPxs7ODm9vb9avXw+kFfoTJ04A8ODBAzw9085nrly5MltzJkfeR3fiNPZxt4izckZXrhIOvsWxsH/Z03iZVLw45vPnYnHtMrJPX/qKr5m4riyBJZfy07osXbJSlDwhPiWeXht70XJVS6zMrdjTdQ8/v/8zZV3Kah3tkXx/J3puK2hTkHIu5ZBScv7OeeJT4v/zerVq1ahSpQpr1qzhxx9/ZPny5VStWpWKFSuyYcMGAAICAmjfvj01atSgcOHC2ZJL6vQknb2MdVgoOmlBrGd5ClTxxsrROlv2/0IeHlguWYjZieNY+1dmTnwfCnVsxfAO14mPf/HmipKXnIs6R61ltVh+bDmf1/uc432O07BEQ61jPUGktU7yB39/f/n4hFLnzp2jQoUKuZ4lWZfMxeiL6Aw6yrqUfWFX35ykj43HEHIZC0Myd22KUqCsO5ZWr/7dIsvvrZToFy9FN/gzklLNGV90Gf13tcfH55UjKYrR++XsL3T5vQt2lnasbruaFqVaaB0JIcQRKeUTF1xUC0Qj1hbWlHUpi4WZBRejLz7REsktqbfuIC6cx2CQ3HUtj3PFotlSPF6JEJj364P1+ZNQwZd5N95nW+XPCPorVdtcipKDpJTM2j+L9uvbU6VIFY71OWYUxeN5VAHR0MMiYm5mTsjdEJJ0uTj4oJSkXgnHMvwqccKRlNK+uJRwMK5utKVK4XR8Nw86D6B/8hysXm/OimlPv0lKUUyZQRoYuGUgw/8eTjvfduzosiNH7tvIbqqAaMzawpqyzmkXxUKiQ0jV58K3bL2e1POXsIy+xR0zVywr+OBYyEh7dFtZ4bTySxK/+YFaZsE0GlWPBQMuko/OvCp5nN6gp8fGHiw6vIhhdYextt3aTA+DpDVVQIyAjaUNZZzLkGJIIfRuKAaDIecOpteTei4Ei/j73LIqhlPlEtjaGVOz4+lse36IxT+7cLOJ4YNF9ZjX8QA5+TYpSm7QG/R029CN745/x/jG45nRcgZmwnQ+lk0naR7nYOVAqYKliE+N59qDa+RI5wadjtSzF7FIiuOGTSlcKxYxqfGozOvWwv7Ev0ingvT5qRkLWv+piohisvQGPZ1/78yqk6uY1HQSAU0CTG4aBFVAjEgh20IUdSxKdGI0t+NvZ+/OdTp0Zy9inpxAhE0Z3Cs4m+TouGZly1D4wn7ue1Tgk7/fZfFbm9TpLMXkSCkZ/NdgVp9azZRmUxjbaKzWkbJEFRAjYG5ujp+fH5UqVeLT7p9ibbAmLCaM2OTYl95X165d+fnnnwHo2bMnZ8+eBb0e3bkQzFISuWFbBo8KBTE3h127drF///6XPkbJkiW5c+fOS2+XXUQRN4qe2U6Ue2V6bXmPFe/9oYqIYlKm/jOVhYcX8lndzxjdcLTWcbJMFRAj8HA499OnT2NlZcXWdVuxsbDh0r1LpOhT0Ol0WdrvsmXL8C1fHt2FUMyT44mwLo1HeadHLY+sFhCjUKgQnme3c9PNj49+b8uaDzdpnUhRMuXbY98yZscYPqz8ITNaztA6zitRBcTINGzYkMuXLhNxMoJubbrR+s3W+Pr6otfrGT58ODVr1qRKlSp8/fXXQFpTeMCAAZQrV44WLVpw+/b/n/pq0qQJB37dgEVCLD8cusw7nZtRvXpVmjdvztWrV1myZAlz587Fz8+PvXv3EhUVRdu2balZsyY1a9Zk3759AERHR9OqVSsqVqxIz549c+b6TBaIQgUpdm4bES5VeHdNezaO2Kd1JEV5rh1XdtD7j960Kt2Kb9t8a1IXzJ/GSPtuauTZ47lnnZ8fzJuXqVV1Oh1btmyhdevWWFtYc/H0RVYHraZWxVosX74cJycnDh8+THJyMvXr16dVq1YcO3aMCxcucPbsWSIjI/H19aV79+4gJTIpCYuEWE7G2DJ66mfs2bMHb29v7t69i7OzM3379sXBwYFhw4YB8MEHHzBkyBAaNGjA9evXee211zh37hwTJkygQYMGjBs3js2bN7N8+fLsfY9egZlzQYqd2kJkmfo0nPEW+0rtpX6fSlrHUpQnXL53mfbr21OucDnWt1+fJ6a/VgXECDwczh3SWiA9evRg//791KpVi0plKxERG8Gff/3J2dNnH13fePDgASEhIezZs4dOnTphbm5O0aJFadasGQCGyNuI1FSihQuxdy/TqFGjR8O+Ozs7PzXH9u3b066ZpIuJiSEuLo49e/bw66+/AvDmm29SqFChp26vFQsPV5wObiOhWj28P3mNcyX2U6F17s+NoCjPEpscyztr3kFKycaOGylgXUDrSNlCFZCMMtlSyG4Zp7TNyN7enhIFSxAfFU9cShzz5s/jjdff+M86f/755xPbyYQERHgYOiyw9nIjJe5WpnIYDAYOHDiAjY1p3MSUkUOlksRv2Ypdq0bEvP0WkWf3UcQnb/wnVUybQRro/Htnzt85z18f/UVp59JaR8o2pn0CLh+wMLOgVKFS1GpUi3kL55Gamnan+sWLF4mPj6dRo0asW7cOvV7PzZs32blzJ/JmJInYIq1tcHAU1KlThz179nDlyhUA7t69C4CjoyOxsf/f06tVq1Z8+eWXj54/LGqNGjVi9erVAGzZsoV79+7lyt/+soq0qMydr9ZTRneOS7U/ICVRr3UkRWH2/tn8fv53ZrWaZfRjW70sVUBMgIOVA31798WzlCd+1dK6+/bp0wedTse7776Lj48Pvr6+dP74Y+pUrIQBuOvsg6VV2k1Jrq6uLF26lPfee4+qVavSoUMHAN5++21+++23RxfRFyxYQHBwMFWqVMHX15clS5YAMH78ePbs2UPFihX59ddfKV68uFZvxQuV6duC492/pN69zeytN0LrOEo+tz9sP6OCRtG2Qls+rf2p1nGynabDuQshWgPzAXNgmZRy2mOvzwWapj+1A9yklAXTX9OTNkshwHUp5TsvOp4xDef+sgzSwNmos+gNeiq6VcTC7LGzj1KiD72M2YN7XLcpRzFfR15m1tmcoOV7+0+1gTQ4vpB/eqygwbKummRQ8rfohGiqfV0NS3NLjvY+ipONk9aRsszohnMXQpiTNl3t64Av0EkI4ZtxHSnlECmln5TSD/gS+DXDy4kPX8tM8TB1ZsIM74LepBpSCY8Jf+J1Q+RtzB/c46aZJ+4+2hcPrdX5dy5HCjWnxvJPuLAum3vWKcoLSCnpuqErkfGR/NTuJ5MuHs+j5cdMLSBUSnlZSpkCrAXaPGf9TqTNoZ5v2VvZ4+7gzp2EO9xPuv9ouYyLQ4SHc5+C2Jd2xzqXJhE0ZhY2FpTYt4b7Zi7YfNyO+Ij7L95IUbLJ4uDFbLq4iZktZ1KjaA2t4+QYLQuIJxCW4Xl4+rInCCFKAN7AjgyLbYQQwUKIA0KI/z3rIEKI3unrBUdFRWVHbk0VdSyKrYUt1+5fQ2/Qg16PIfQKKViS6F4SJyfTGowtJxWu4MrNeesomnqNC/W6ocY7UXLDhTsXGLZtGK3LtGZgrYFax8lRpnKioyPws5QyY7eaEunn5D4A5gkhnto3Tkq5VErpL6X0d3V1zY2sOcpMmFGiYAlSDanciL2B7sp1zHTJRNqVwt1T9cp+XPWB9dnecgbVr//O0Y9max1HyeN0Bh2df++MraUty99ZbnKj674sLQtIBFAsw3Ov9GVP05HHTl9JKSPSf14GdgHVsj+icXKwcsDVzpXI+EiS46OJFB64lzGy2QSNSMvNg9nl0pYqq0cS/pOJjv2lmITAvYEcijjE4jcXU9SxqNZxcpyWBeQw4COE8BZCWJFWJDY+vpIQojxQCPg3w7JCQgjr9N8LA/WBs49vm5d52rhhqYcrBc2w8vbAyvRHRcgxFpaCMnu+JcKsGOLjj0iNjtE6kpIHHY44zMTdE/mg8ge8X/F9rePkCs0KiJRSBwwAtgLngJ+klGeEEBOFEBl7VXUE1sr/9jeuAAQLIU4AO4FpUkqTLSAPh3OvWLEiVatWZfbs2c+flVBKuHwdzxhBkqWBVKu0azv16tXL0vGvXr366EZBgODgYAYNGpSlfRkrL98CXJn4A+4p1zjbIm/9bYr2knXJdPm9C+4O7ix8faHWcXKPlDLfPGrUqCEfd/bs2SeW5TZ7e/tHv0dGRsrmzZvLcePGPXN9/c1IKQ8flmHHbssLURfkwesHZbIuOcvH37lzp3zzzTezvP2zGMN7+7hfK46VEuTl6T9pHUXJQ8btGCcJQP558U+to+QIIFg+5TPVVC6i5xtubm4sXbqUhQsXIqV8chj3hQshIpxNR87xfr/3GNptKO2btCfsQRgODg4AdOzYkc2bNz/a58NJpq5evUrDhg2pXr061atXfzQXyMiRI9m7dy9+fn7MnTuXXbt28dZbb2EwGChZsiT37/9/F1gfHx8iIyOfOfS7sWscNI5jljVxHt2H5EtP3k+jKC/r9O3TTP1nKh9W/pDXfV7XOk6uUt12Mhj812CO38rem8783P2Y1/rlBmksVaoUer2e27dvs2HDhv8fxj0piXr+/jQvVoJkpyKcOHGU06dPY13YmhuxNx5t36FDB3766SfefPNNUlJSCAoKYvHixUgp+fvvv7GxsSEkJIROnToRHBzMtGnTmDVrFps2pU3KtGvXLgDMzMxo06YNv/32G926dePgwYOUKFGCIkWKPHPod2PnXMSSU4t+wKJ3NcKad6XMlb9RvQ+UrNIb9PTc2BMnG6eX/n+eF6gWiJHbtm0b33//PX5+ftT29+fuvXscvpFIITcratWqhbe3N+727liZW2GQBqSUvP766+zcuZPk5GS2bNlCo0aNsLW1JTU1lV69elG5cmXat2//n6Hbn6VDhw6sW7cOgLVr1z4aR2v79u0MGDAAPz8/3nnnnUdDv5uCxr3K8mv9OZS5FsSlkUu1jqOYsC8PfcnBiIPMbz2fwnaFtY6T61QLJANj+QZx+fJlzM3NcXNzQ0rJl19+yWtNm6I/dYYEaYsoV47gI7uxt7cH0loKxQqk9YiOSojCzd6NJk2asHXrVtatW0fHjh0BmDt3LkWKFOHEiRMYDIZMDdtet25dQkNDiYqK4vfff2fs2LGAaQ/9DtBmc2/2FfmJqrOGk9zzdax9jHeASMU4Xb1/lTE7xvB6mdfpVKmT1nE0oVogRiYqKoq+ffsyYMAAhBC89tprLF68mIQLl0BKTjxIRZglPLFdQZuCCCGIiIkgVZ9Khw4dWLFiBXv37qV169ZA2iRUHh4emJmZsWrVKvT6tPsyHx/WPSMhBO+++y5Dhw6lQoUKuLi4AM8e+t1UFHASpCz6BmHQc/2NPuoudeWlDdwyEIFgyVtL8vwNg8+iCogReDgjYcWKFWnRogWtWrVi/PjxAPTs2ZMKpUpR+923qdjxA8ZOGoROp3tiH0IIBAKDNHAj9gatWrVi9+7dtGjRAqv0m0T69evHypUrqVq1KufPn3/UgqlSpQrm5uZUrVqVuXPnPrHvDh068MMPPzw6fQU8c+h3U9K0RynWV5+GT+hf3Jj2vdZxFBOy8cJGNl3cRECTAIo75d/Wq6bDuec2kxzO3WBAd+IMqXpBqo8vBZyeX/PDHoQRGR9JhcIVsLeyz6WQT2f07y0QedPA5eKNqSRPY3/tLGaeHlpHUoxcQmoCvot8cbBy4FifY1iaW2odKccZ3XDuSuakht3EQp/MPccSLyweAB6OHliYWRAWE0Z++nKQVUU8zLgxcTkW+iSuvJ23B75Tskfg3kCuPbjGojcW5Yvi8TyqgBgxmZSEedQt7gpnXEs5ZmobCzMLPB09iUuJ+8+Q78qzvTeyLKtLfUHpY79w98ctWsdRjNjF6IvM3D+Tj6p8ROOSjbWOozlVQMA4v6lLSeql6xgww1DUC8uX+KJT2K4wtha2hMeEY5DPGRIlBxnle/oMQkDjP4ZxnvKk9h0AiYlaR1KMkJSSgVsGYmNhw8yWM7WOYxTyfQGxsbEhOjra6D7w9NH3sEqM4Y5VUVzcX26kRCEEXgW8SNYnExWf+3OgSCmJjo42qS6+ZXyt2P/RVxSJu8y1PoFax1GM0C/nfmHbpW1MbjoZdwd3reMYhXx/ET01NZXw8HCSkpI0SvUUBgOG8BvopDm4u2NlnbUugpFxkaToU/As4ImZyN3vCjY2Nnh5eWH5Mk0njSUkwF9uH/N2wjo4cRLLyuW1jqQYicTURMotLIezrTPBvYOxMMtft9A96yJ6/noXnsLS0hJvb2+tY/xHRK8APJdN4Mv39zJwne+LN3gGXaQOv6/9+LT2p8x5bU42Jsyb7OzAduEs4rptIu69fhS7GKSGOVEAmLV/FmExYfzw3g/5rng8T74/hWVsZFg4Lt/OYIP1+3T5psEr7atykcp09+vOwkMLCb0bmk0J87bWXYrwg+9UioXu5N5Xq1+8gZLnRcREMG3fNNr5tqNRiUZaxzEqqoAYmSudRoHBQMqk6RQo8Or7m9h0IlbmVozcPvLVd5YPCAGtf+3NIVELMewziFGTT+V3o4JGoTfomdFihtZRjI6mBUQI0VoIcUEIESqEeOITTgjRVQgRJYQ4nv7omeG1LkKIkPRHl9xNnjMSdh2i1L4fWOs+hLaflcyWfXo4evB5/c/55dwv/HP9n2zZZ17nU86MI10XUjApkrBP1AX1/Oxg+EFWnVzF0LpD8S5kXKe6jYFmF9GFEObARaAlEE7aFLedZIaZBYUQXQF/KeWAx7Z1BoIBf0ACR4AaUsp7zzvm0y6iGw0puVq8ITbhoVz/+yK1WmRD8yNdfEo8ZReWxauAF//2+DfXL6ibooQE2OzalTaJa7A4fwazsmW0jqTkMikl9b6tx9X7V7k44CKO1pm7FysvMsY70WsBoVLKy1LKFGAt0CaT274G/C2lvJteNP4GWudQzlwRuXA9JcP3san25GwtHgD2VvZMaTaFQxGHWH9mfbbuO6+yswPzGVNJllaEdxymdRxFA2tOr+FA+AECmwXm6+LxPFoWEE8gLMPz8PRlj2srhDgphPhZCFHsJbdFCNFbCBEshAiOisr9eyIyJSkJMfJzTppV5fWfuuXIIT6u8jGV3CoxdudYUvWpOXKMvObdfh58X2wMxY9tIPGP7VrHUXJRQmoCI7aPoLpHdbr45Ykz5DnC2M9l/AGUlFJWIa2VsfJldyClXCql9JdS+ru6umZ7wOwQMuhL3BKucarrHDyLm+fIMczNzJnafCqhd0NZdnRZjhwjrxECanw/mEuUIqbHYHjKKMhK3jRr/yzCY8KZ99o8dcr3ObR8ZyKAYhmee6Uve0RKGS2lTE5/ugyokdltTYUu6h5u305ll+3rtF3ULEeP9abPmzQo3oCJeyYSnxKfo8fKK+o0seH3BrMpEnWGe9O+1jqOkgtuxd1ixr4ZtK3QloYlGmodx6hpWUAOAz5CCG8hhBXQEdiYcQUhRMaxtd8BHk66vRVoJYQoJIQoBLRKX2ZyznSejqP+PrpJU8npkT+EEExvMZ1bcbeYd8A4Zl80Be1WtWGHWXMsJ30Bd+9qHUfJYRN3TyRZn8zU5lO1jmL0NCsgUkodMIC0D/5zwE9SyjNCiIlCiHfSVxskhDgjhDgBDAK6pm97F5hEWhE6DExMX2ZSEkIiKPvXfLYV/pDmQ6vmyjHrFatHm3JtmL5vOncS7uTKMU1diZKC093nYpdyn5uD1IdKXnbhzgWWHllK3xp98XHx0TqO0cv3Y2Fp6ah/LyodWcmp9Reo0S73+pifjTpL5cWV1RAnLyE2FjYV6U7bpB+xvHQB4V1S40RKTnhv3Xtsv7yd0EGhuNm7aR3HaBhjN9587c4/56l65Fu2le6Xq8UDwNfVl65Vu7Lo8CKu3b+Wq8c2VY6OwISJ6KUZ17t8oXUcJQfsu76P387/xuf1P1fFI5NUAdFIWOcxxGNPhR/GaHL8gCYBCATjd43X5Pim6P2hXvxQeAgl9v5A6qFjWsdRspGUkuF/D8fDwYMhdYZoHcdkqAKigWvrDlDtyq/sqTmM0nW06VpczKkYA2sN5PsT33Mq8pQmGUyNuTkUWziCO7gQ2Xk45KPTv3ndb+d/49/wf5nYdCL2VvZaxzEZqoDkNimJ7T+S28KN2uuGahplVMNRFLAuwOgdozXNYUpee9+J1aW+wOtCEIkbtmkdR8kGqfpURm4fmXZq16+r1nFMiioguez03L+pFL2bo298gau3g6ZZnG2dGdlgJJsubmLvtb2aZjEVQkCtFZ+k3VzY93PQ67WOpLyiZUeXEXI3hOktpqu5Pl6SKiC5SBokcnwAEebFaLSql9ZxABhUexBFHYsyMmik0U3ra6zqNLLi95qBFIk8ScxXP2gdR3kFscmxBOwOoFGJRrzp86bWcUyOKiC56J/xf1M57l+udhqNXSFrreMAYGdpx/jG49kftp/NIZu1jmMy3vyuPYeoiX70WDCm6ZCVlzJr/yxux99mZsuZCDX75EtTBSSX6FIl9rMncNOiGHW+zpkBE7Oqm183ShcqzZgdYzBIg9ZxTEJ5XzP2vDGdQnHhRE9ZonUcJQtux99m9r+zaefbjlqetbSOY5JUAcklQaO3Uz1xPze7jXJOK/gAACAASURBVMbczjhaHw9ZmlsyselETkae5KczP2kdx2R0WtqUHWbNsZgZCHFxWsdRXtK0f6aRqEtkUtNJWkcxWaqA5ILkJInLlwHcsixGtQXG1fp4qGOljlR2q8wXO79Qw71nkqcnnP9oCk7JUdwcOV/rOMpLCI8J56vDX9G5amfKFy6vdRyTpQpILvhr2Hb8k/dzt88ohI1xtT4eMhNmTG42mdC7oaw88dKj5udbH8yvzZ+W7+D49Uy499wJMRUjMmn3JAzSwPjG6kbaV6EKSA6Lj5N4LA3gtpUXFWZ21zrOc71d9m1qe9Zmwu4JJOnUheHMKFgQbg+YhJ0uhrBBM7WOo2RC6N1Qlh9bTp8afShZsKTWcUyaKiA5bNPg7dRK3U/MgNFG2/p4SAhBYPNAwmPCWRKsLgxn1vuTq7DBugOFV89H3orUOo7yAgG7ArAyt2JMI22GEcpLVAHJQQ/uS0p8N4EoGy/KBBp36+OhZt7NaO7dnMC9gcQmx2odxyTY2UHc8AlYGpIJ+yRQ6zjKc5yKPMXqU6sZVHsQ7g7uWscxeaqA5KDfBgRRR7+PxMGjwdq4Wx8ZTWk2haiEKOYfVBeGM+v9sWVZb98V9w1LkNeuax1HeYYvdn6Bo7Ujn9f/XOsoeYKmBUQI0VoIcUEIESqEGPmU14cKIc4KIU4KIYKEECUyvKYXQhxPf2x8fFutRd2WlFsTwB1bL4oHmEbr46HaXrVpU64NM/fP5G6iyc3TpQlra5BjxyElhPVW3UKN0aGIQ2y4sIFhdYfhbOusdZw8QbMCIoQwBxYBrwO+QCchhO9jqx0D/KWUVYCfgRkZXkuUUvqlP97ByPzSL4i6hn3oho0yqdbHQ5ObTSY2OZYZ+2a8eGUFgPafFWdNgb4U3bYCw4UQreMojxmzYwyF7QozuM5graPkGVq2QGoBoVLKy1LKFGAt0CbjClLKnVLKhPSnBwCvXM6YJeFhkiq/BnDXzhP3MT20jpMlldwq8WGVD1lwcAE3Y29qHcckWFqC/ZTRJGNNWK8JWsdRMth5ZSfbL29ndIPROFo7ah0nz9CygHgCYRmeh6cve5YewJYMz22EEMFCiANCiP89ayMhRO/09YKjoqJeLXEmrf9kB/XkPuRI07r28biAxgGkGlKZsneK1lFMxnufFGGtc3+89q5Bf/aC1nEU0iaLGrNjDJ6OnnxS8xOt4+QpJnERXQjxEeAPZOxoXyJ9jt4PgHlCiNJP21ZKuVRK6S+l9Hd1zfnJmy6FSmpuDuCevScun5tm6+Oh0s6l6VmtJ0uPLOXKvStaxzEJ5ubgOn0YSdhwvZe6FmIMNods5t/wfxnXeBw2FjZax8lTtCwgEUCxDM+90pf9hxCiBTAGeEdKmfxwuZQyIv3nZWAXUC0nw2bWrwN20IB/EKNNu/Xx0NhGYzE3M2fCbnVKJrPe6u7Getf+FN+/htRT57WOk68ZpIGxO8ZSulBpuvkZ5zBCpkzLAnIY8BFCeAshrICOwH96UwkhqgFfk1Y8bmdYXkgIYZ3+e2GgPnA215I/Q2iIpO7WAO47eFLwM9NufTzkWcCTATUHsOrkKs5Gaf4WmwQzM/CYndYKudZrstZx8rX1Z9ZzIvIEE5tOxNLcUus4eY5mBURKqQMGAFuBc8BPUsozQoiJQoiHvapmAg7A+se661YAgoUQJ4CdwDQppeafbr88bH2MMs2eV88yosEI7C3tGbdznNZRTEarj9z41b0/3gfXkHJStUK0oDPoGLdrHJXcKtGxUket4+RJIj/NQufv7y+Dg4NzZN+hIZKbZRtT2eEyBaNCwSZvnWudsGsCAbsDONzrMP5F/bWOYxJ2/XSbmh28Cfd/l3KH1cyFue3bY9/SY2MPfu/wO23Kt3nxBsozCSGOpF9z/g+TuIhuCn7uv5OG7MVs9Kg8VzwAhtQdgoutC2N3jNU6islo3N6NDV4DKBO8huQTqhWSm5J1yQTsCqCWZy3eKWd0t4nlGaqAZIOQi5J6f6dd+ygwJG9c+3hcAesCjGowiq2XtrL76m6t45gEIaDY/LRrIVd7qh5ZuenrI18TFhPGlGZT1FS1OUgVkGzwc/+dNMrDrY+H+tXsR1HHoozZMYb8dOrzVTR415UNXgPwUa2QXBOfEs+UvVNoWrIpzb2bax0nT1MF5BVdvCCpvz2ABw5F82zr4yFbS1vGNRrHvrB9bAnd8uINlEetkERsVSsklyw4uIDb8bdV6yMXqALyin4ZkD9aHw91r9adUoVKMTpoNAZp0DqOSVCtkNxzP+k+M/bP4K2yb1G3WF2t4+R5qoC8ggvnJfW2T+CBQ1Ech/TUOk6usDS3ZGKTiZyIPMH6M+u1jmMShIDiC4aRgJ1qheSwWftncT/pPpObqvtvcoMqIK/g5wG7aMyefNP6eKhjpY5UcqvEFzu/QGfQaR3HJNT/nysbiqW3Qo6f0zpOnhQZF8m8A/PoWKkjVd2rah0nX1AFJIvOn5M0CArIV62Ph8zNzJncdDIhd0NYeXyl1nFMghBQYv5n6a0Q9e04J0z9ZypJuiQmNFHD7uQWVUCy6JeB6a2PMfmr9fHQO+XeoZZnLSbsnkCSLknrOCah/v9c+cOrH2WOrCX5pBqpNztdf3CdxcGL6erXlbIuZbWOk2+oApIF587x/62Pwfmr9fGQEILAZoGExYSx+PBireOYhIz3hVzpqYbIz06TdqddWxrXWA23k5tUAcmCh9c+zEePzJetj4eal2pOi1ItCPwnkJjkGK3jmIT677rxh9cn+Bz+keTTatbC7BASHcKK4yvoW6MvxZ2Kax0nX1EF5CWdPQsNdwQQ4+CBw5BeWsfRXGCzQO4k3GHOv3O0jmIShIBi84aRghWXe6hWSHYYv2s81hbWjG44Wuso+Y4qIC/p5wG7aMLufNfz6llqetakbYW2zP53NlHxuTPjo6mr39adTZ598Tn0A8lnL2kdx6SduHWCNafXMLj2YIo4FNE6Tr6jCshLOHMGGu1UrY/HTW42mYTUBAL3BmodxWR4zv+cVCy51F21Ql7FFzu/wMnaiWH1hmkdJV/KVAERQrgJId4VQvQXQnQXQtQSQuS74vOw9ZHfr308rnzh8nSt2pWvgr/i2v1rWscxCfXaerDZszdlD35P8nk1XXBWHAg/wB8X/+Dz+p9TyLaQ1nHypecWASFEUyHEVmAz8DrgAfgCY4FTQogJQogCOR9Te6dPQ6NdE4hx8MB+sGp9PG58k/EIhJr69iUUnTcCHRaEdlctt6wYs2MMbvZuDKo9SOso+daLWhFvAL2klDWllL2llGOllMOklO8AVYFjQMusHlwI0VoIcUEIESqEGPmU162FEOvSXz8ohCiZ4bVR6csvCCFey2qGzPp5wC6asiut9WFrm9OHMznFnYrTr2Y/Vp5Yqaa+zaR67YryZ9FelP33O5IvXNU6jkkJuhzEjis7GN1gNA5WDlrHybeeW0CklMOllNef8ZpOSvm7lPKXrBxYCGEOLCKtZeMLdBJC+D62Wg/gnpSyDDAXmJ6+rS9pc6hXBFoDX6XvL0ecPg2Nd6vWx4uMajAKe0t7vtj5hdZRTIbHvBEYMCOk+1Sto5gMKSVjdoyhWIFi9PHvo3Ucoyd1eu4dypku45m9BrJKCOGU4XlJIUTQKx67FhAqpbwspUwB1gKPzzvZBng4VsbPQHORNj5zG2CtlDJZSnkFCE3fX45Y31+1PjLD1d6VYfWG8eu5XzkUcUjrOCahbnsvthTtQdn9K0i6+NTvaspj/rj4BwcjDjK+8XhsLNS1yBc5NmItBWqX59zK7P8/mdkL4f8AB4UQbwghegHbgHmveGxPICzD8/D0ZU9dR0qpAx4ALpncFgAhRG8hRLAQIjgqKmvdTPvfnUi8k2p9ZMaQOkNwtXNlVNAoraOYjCJz0s7eXlStkBcySANjd4zFx9mHLn5dtI5j9KROj/OiiYRYVcSn0xNTmr+yTBUQKeXXQE9gAzARaCSl/CPb0+QAKeVSKaW/lNLf1dU1S/twWz0f+9XLVOsjExytHRnTcAw7ruxg++XtWscxCXU7FOcvj+6U37ecpJCwF2+Qj607vY5Tt08xockELMwstI5j9I6NWEvJ5Ivc7DUeC6vs7zib2VNYHwPfAp2B74A/hRCvOl5yBFAsw3Ov9GVPXUcIYQE4AdGZ3Db7VK4Mb7yRY7vPa/r6pw0pMSpolJr6NpPc5qa12C50n65xEuOVqk/li51fUKVIFTpU6qB1HKP3sPVx3qoyDee8myPHyGxJags0kFKukVKOAvqSVkhexWHARwjhLYSwIu2i+MbH1tkIPGyntgN2yLRPpI1Ax/ReWt6AD6BOuhsJawtrJjSZQPCNYH4996vWcUxCnQ4l2ObRlfL/fENiaM59FzJlK46v4NK9S0xpNgWz/Hcb2kvL6dYHgMjqN0QhhFX6xe+sH1yIN0i7lmIOfCulnCKEmAgESyk3CiFsgFVANeAu0FFKeTl92zFAd0AHDJZSvnCSbn9/fxkcHPwqkZVM0hv0VF5cGYM0cLrfaXW6IRMOrbtCtY5lOd3gE6rtXaB1HKOSmJqIz5c+FHcqzr7u+9Rc5y8gdXquOfiSJK0pE3v8lQuIEOKIlPKJiygvupFwrBDC+akBpUwRQjQTQryV1VBSyj+llGWllKWllFPSl42TUm5M/z1JStleSllGSlnrYfFIf21K+nblMlM8lNxlbmbOlGZTuBB9QU06lUm1Oniz3aMzFf5ZSuKlG1rHMSqLgxcTERvB1OZTVfHIhNxofcALWiBCiDbA50AScBSIAmxIO2XkB2wHAqWUJjGKnmqB5C4pJXWW1+FG7A1CBoaoLpeZcPiny1TrUJaTDQdQfc+rdnTMG2KSYyg1vxQ1itZg60dbtY5j9NJaHxVJklbZ0vqALLZAgHZSyvrAVuAMaaeaYoAfgFpSyiGmUjyU3CeEYGrzqYTHhPPV4a+0jmMSar5fiiCPj6mw92sSLt/SOo5RmHdgHtGJ0UxppgaezIxjI9dRMvlCjrc+4MUtkLNAC2AL0PTx16WUd3MuWvZTLRBttFrViqM3j3L508sUsM4XQ6e9kiPrQvHrWI5jjQbjv3u21nE0dSfhDqXml6Jl6Zb88n6WBr3IV3Ki9QFZb4EsAYKA8kBwhseR9J+K8kKBzQOJToxm9v78/WGYWTU6lGGHx4f47llM/OVIreNoavo/04lLiWNS00laRzEJudn6gBePhbVASlmBtB5SpTI8vKWUpXI8nZIn+Bf1p51vO2b/O5tbceq0TGa4zBmLNcmc6TZL6yiaiYiJYOHhhXxc9WN8XR8fJk95nNTpcV6Ys/d9PC6zd6J/ktNBlLxtSrMpJOmSmLh7otZRTEL1jmXZ5dGJinu+Iv7Kba3jaGLynsnoDXoCGgdoHcUk5HbrA9SMhEouKetSlj41+rD0yFIu3LmgdRyT4DJnLLYkcqp7/ptv/tLdSyw7tozeNXrjXchb6zhGT4vWB6gCouSicY3HYWtpy+gdo7WOYhL8OpZnt0dHKu9aSNzVO1rHyVXjd43H0sySMQ3HaB3FJGjR+gBVQJRcVMShCJ/X+5xfz/3K/rD9WscxCWmtkAROdcs/rZBTkadYfWo1g2oPwsPRQ+s4Rk+r1geoAqLksqF1h+Lu4M7wv4ergRYzoUpHX/5xb0/lXV8SezVa6zi54oudX1DAugCf1/9c6ygm4ejInzRpfYAqIEous7eyZ0KTCewP28+GCxu0jmMSCs35AgfiONltrtZRctzB8INsuLCB4fWG42z71FGUlAy0bH2AKiCKBrpX6075wuUZsX0EqfpUreMYvcqdKvGPezuq7FpA7DWTunf3pY3ZMQZXO1c+rfOp1lFMQvDnP+GdfF6T1geoAqJowMLMguktpnMx+iLLjy3XOo5JKDhnHI7EcrzbfK2j5Jigy0EEXQliTMMxOFg5aB3H6BlS9bh8NZELGrU+QBUQRSNvl32bBsUbELArgLiUOK3jGL1KnSqz3/09quyaT8z1+1rHyXZSSsbsGEOxAsXo499H6zgm4eDgNZRKPk9Uv3GatD5AFRBFI0IIZracSWR8pBriJJOcZn2Bk3yQJ1shv53/jYMRBxnfeLwatTkTdImpeHwzgXM2ftSb+Z5mOTQpIEIIZyHE30KIkPSfhZ6yjp8Q4l8hxBkhxEkhRIcMr30nhLgihDie/vDL3b9AyQ51vOrQzrcdM/fPVEOcZELFD/04UKQNVXbO48H1B1rHyTY6g45RQaPwdfWli1+XF2+gcGjA95RMDeXe4ImYWWjXDtDqyCOBICmlD2mDNY58yjoJQGcpZUWgNTBPCFEww+vDpZR+6Y/jOR9ZyQmBzQJJ1ierIU4yyWn2OArK+xzrlndmLPz22LdcjL7I1OZT1cyVmZAal0yJ7ydyyq4WdadkeT6/bKFVAWkDPJymbiXwv8dXkFJelFKGpP9+A7gNuOZaQiVX+Lj4qCFOXkKFD6tzqMjbVN05l3vXYrSO88riU+IJ2BVA/WL1ebvs21rHMQkH+yzHU3edhJGTEGbazs6oVQEpIqW8mf77LaDI81YWQtQCrIBLGRZPST+1NVcIYf2cbXsLIYKFEMFRUWruK2P0cIiTUUGjtI5iEgrNG08heY/DXRZqHeWVzT84n5txN5nWYpqaqjYTku8nUmbdFI45NqTWmJZax8m5AiKE2C6EOP2UR5uM68m025GfeUuyEMIDWAV0k1Ia0hePIm2OkpqAMzDiWdtLKZdKKf2llP6urqoBY4zc7N0YUX8Ev53/jX3X92kdx+j5dKzB0aJvUmP3bG5fitU6TpZFJ0Qzfd903in3Dg2KN9A6jkk41H0J7vob6Mdr3/qAHCwgUsoWUspKT3lsACLTC8PDAvHU8aqFEAWAzcAYKeWBDPu+KdMkAyuAWjn1dyi5Y0idIRR1LMpn2z7D8Oh7gvIsheePw4W7HOy8SOsoWRa4N5C4lDgCmwVqHcUkJNyOo/yGaQQXbEGNoY21jgNodwprI/Cwu0UX4IkxLYQQVsBvwPdSyp8fe+1h8RGkXT85naNplRxnb2XPlGZTOBhxkLWn12odx+gVb1eLU16tqbt/FmHnTO8+mmv3r7Hw8EK6Vu1KRbeKWscxCUe6LcTVcBvzwEkYy9k+rQrINKClECKEtDnXpwEIIfyFEMvS13kfaAR0fUp33R+FEKeAU0BhYHLuxldyQueqnanuUZ2R20eSkJqgdRyj57ZwPIWJ5mCXr7SO8tLG7RqHmTAjoEmA1lFMQmz4AyptmcFBlzeo9kkdreM8okkBkVJGSymbSyl90k913U1fHiyl7Jn++w9SSssMXXUfddeVUjaTUlZOPyX2kZTS9L6CKU8wE2bMfW0uYTFhzPk3/wxfnlVF2tThfPFWNDo8i5Dj8VrHybSTkSdZdWIVg2oNophTMa3jmIRjXeZRSN7DbpZxdXdXd6IrRqVRiUa8V+E9pv0zjRuxN7SOY/TcvhqPG1Ec6mo6rZBRQaNwsnFiZIOn3f6lPO5u6F2q7pjDfvd3qdy1htZx/kMVEMXozGgxg1RDKmN3jNU6itFzfrMeF0u9RusT0zj1j/Hfnb7r6i7+DPmT0Q1GU8j2iQEolKc49tEsHIml8MIJWkd5giogitEp7VyaQbUG8d3x7zh686jWcYyex7eBuHCXM92Ne0wxKSUjto/Aq4AXA2oN0DqOSYgIvkmdg/M5VLIDZdtW1jrOE1QBUYzS2EZjcbFzYejWoWrmwhdwbFydc5Xf562QOQRvjtQ6zjP9fPZnDkUcIqBxALaWtlrHMQkhH0/EihS8VkzSOspTqQKiGCUnGycmNpnI7mu7+f3871rHMXolvp+EDUlc6xOIMdbbZF0yI7aPoJJbJbr6ddU6jkkI3RJCg/PfcKBKb7yalNE6zlOpAqIYrV41euHr6svwv4eTrEvWOo5Rs/Mry4V63Xk7YjG7V17VOs4TFh5ayJX7V5jVchbmZuZaxzEJkb3Gkow1FX78Qusoz6QKiGK0LMwsmN1qNpfuXWLhIdMf9ymn+awahxRmPBgSgF6vdZr/F50QzeS9k3mt9Gu8VuY1reOYhFMrj1A/4ieONBpK4UruWsd5JlVAFKPWukxrWpdpzaQ9k4iKV4NhPo9VKS+uvjWQt+9/z8apZ7SO88ikPZOISY5hVqtZWkcxCVJC8pCRRAsXqq0ernWc51IFRDF6s1vNJi4ljoBdAVpHMXplvx1JgrkjtlPGkmAEN/OHRIew6PAielTrQSW3SlrHMQnB07bjf287Z/43FkfPAlrHeS5VQBSj5+vqS1//viw5soQTt05oHceoicIu3Ok6nNZJv7P+swMv3iCHjdg+AhsLGyY2Na47qI2VQWfAftJIwi1KUGflJ1rHeSFVQBSTMLHpRArZFGLgloGqW+8LlJw3mPtWbnh/M4qo29q9V3uu7eG3878xov4I3B2M9zy+Mfln8M/4Jh7heo+JWDk+c5ojo6EKiGISnG2dCWweyN7re9VovS/i4EDSsLE00u/i155/ahLBIA18tu0zPB09GVp3qCYZTE3Cg1SKLRlDqG0l6nz5odZxMkUVEMVk9KjWg+oe1Rn29zDiUtT4mc/jPr4PkU4+NPxjOKHndbl+/DWn1hB8I5gpzaZgZ2mX68c3Rbs/Woq3PpTk8VMxszSNrs6qgCgmw9zMnIWvL+RG7A2m7JmidRzjZmWF5ZwZ+HKOnR9+k6uHTkxNZFTQKKq5V+Pjqh/n6rFN1a3z96m1aTynXJtS8fM3tY6TaaqAKCalbrG6dK7amdn/ziYkOkTrOEbNuVsbrhZvRJuj4zkcFJNrx53z7xzCYsKY3Wo2ZkJ9xGTGifcnU4i7OC2bg9HMFpUJmvzrCiGchRB/CyFC0n8+dVhOIYQ+w2RSGzMs9xZCHBRChAoh1qXPXqjkE9NbTMfGwobBWwdrHcW4CYHrqjm4EcW5zlNzZYiT8JhwAv8J5N3y79LUu2nOHzAPOL/5Ek1PLeBwxW4Uf8fvxRsYEa2+HowEgqSUPkBQ+vOnScwwmdQ7GZZPB+ZKKcsA94AeORtXMSbuDu6MbzyeP0P+ZNPFTVrHMWr2jWoQWvdj3r8xlw3zr+b48UZsH4HeoGd2K+MeGdhYSAlR3T4nBSvKrTe9iVW1KiBtgJXpv68kbV7zTEmfB70Z8HCe9JfaXskbBtYeSPnC5Rn812CSdElaxzFqpdZMASGQo0cTl4N9D/Zd38fqU6sZVm8Y3oW8c+5AeciBmXtoGPUrp94YScEKHlrHeWlaFZAiUsqb6b/fAoo8Yz0bIUSwEOKAEOJhkXAB7kspH3YtCQc8n3UgIUTv9H0ER0WpoTDyCitzKxa0XsCle5fU9LcvYFaiGHc6f8a7iWtYNfBQjhxDb9Az6K9BeDp6MqrBqBw5Rl6TmmzAcfxQblp4UeNH0+zqnGMFRAixXQhx+imPNhnXk2l3hT3r7GwJKaU/8AEwTwhR+mVzSCmXSin9pZT+rq6uL/+HKEarZemWvFv+XSbvmczV+1e1jmPUvL4cwX2bIvitHMKVy9l/MWTF8RUcvXmUGS1nYG9ln+37z4uCuv1ApaQj3Bg4DauCptnVOccKiJSyhZSy0lMeG4BIIYQHQPrP28/YR0T6z8vALqAaEA0UFEJYpK/mBUTk1N+hGLd5redhJszUHeov4uiIYeJk6sr9bOy4Olt3fT/pPqODRlO/WH06VeqUrfvOqyIvx1N57WguONWk+kzTfc+0OoW1EeiS/nsXYMPjKwghCgkhrNN/LwzUB86mt1h2Au2et72SPxR3Ks6EJhPYdHGTmnjqBZw/606EZ03ePzyMPZuyr1vvxN0TuZNwhwWvL0CYUBdULR36XyCeMgK7JXMR5qbb1Vmr5NOAlkKIEKBF+nOEEP5CiGXp61QAgoUQJ0grGNOklGfTXxsBDBVChJJ2TWR5rqZXjMqndT6lapGqDNwykNjkWK3jGC8zM1zWLKIIkVzuNgldNtygfi7qHF8e+pKe1XtS3aP6q+8wHzi69iKvnZrJ0YofU6xjfa3jvBKRn5r9/v7+Mjg4WOsYSg44EH6Aesvr8WntT5nbeq7WcYzalRa98Ar6jvVjT/LBpApZ3o+UktY/tuZg+EFCBobgaq+uMb6IXic55NyainEHMA+5gH1p0xhkUghxJP169H+YbttJUTKo41WHPjX6sODQAo7dPKZ1HKNWcnUgSRYOFJ06kFs3s/4F8uezP7Pt0jYmNJmgikcmBQ34jbqx2wj9eKLJFI/nUQVEyTMCmwdS2K4wfTb1QW8wojldjYxwcyVx9GSa6INY2/6XLO0jJjmGwVsH4+fuR/9a/bM5Yd4UHZZAhW+GcMm+MtWW5Y33TBUQJc8oZFuIua/N5fCNw3x95Gut4xg1t3F9uenux3v7hrLjj/iX3n78zvHcjL3JkjeXYGFm8eINFPa/FUgxw3VYuAhhmTfeM1VAlDylU6VOtCjVglFBo7gZe/PFG+RX5uY4/7iQ4oRxpfN4kl7iZv5jN4+x4NAC+tToQ22v2jmXMQ85sCqEVidncqziR5Tu2lDrONlGFRAlTxFC8NUbX5GsS6b/n/3VvSHPYd2sPuFv9qHr/bms/PRoprYxSAOfbP6EwnaFCWwemMMJ84akRImhzyekCBvKb5yhdZxspQqIkuf4uPgwockEfjv/Gz+f/fnFG+RjXj9MI8bGjZrf9CLk3Iv79X5z5BsORhxkVstZFLJ96iDaymP+7LiSeolBhA+cjm0p0xvv6nlUAVHypM/qfUYNjxr0/7M/dxLuaB3HeBUsCAu+pLo8yva35z93yPfb8bcZGTSSJiWb8FGVj3Ivowm7sCeSJhuHcsG1ARXm9tY6TrZTBUTJkyzMLFjRZgX3k+4z+C81b8jzFOrZlmtV3qbzpXGsnXrlmesN2zaM+JR4vnrjK3XHeSYYDBDWdjD2xOP621Iwm2wykAAAFD1JREFUy3sft3nvL1KUdJWLVGZ0w9H8eOpHNW/I8whB8T8WIczNcBvXl7DrTzZDtoRsYdXJVYyoP4IKrlm/+TA/+WvQn7S4s5Yz/xuDc/28+Z6pO9GVPC1Fn0KNpTW4m3iXM/3OUNCmoNaRjNadCYsoHDCA+ZWWMuhkr0czq8Ykx1Dpq0o4WDlwrM8xrC2stQ1qAq6eisWsaiX0do6UjD6KsDbtSVPVnehKvmRlbsWKNiu4FXeL4duGax3HqBX+4hPCyjan++mh/Dbn/09ljdw+kvCYcJa/s1wVj0wwGOD0a0PxlOHYrvrG5IvH86gCouR5/kX9GV5vOMuOLWNr6Fat4xgvMzM8t36LmbnAbUQ3boQb2HNtD4uDF/Np7U+pW6yu1glNwuZ+m3nr5jJOv/457u/m7fdMncJS8oUkXRI1ltbgXuI9Tn1yChc7F60jGa1bU1fgPro78yvOYFGfb9AZdJz65JSaKCoTLh2Kxr52JRIc3fC+fQhhkzdabOoUlpKv2VjY8ON7P3In4Q59N/dVNxg+h/vIrlyt/BbXio4i5G4I37z9jSoemaDXw9U3+uFM9P+1d+fRUVXZHse/OwMzyDwPSSsNImOIKAi0BBHQhoBEjIog4GO1Qmtro8BrXaKNPeDc2O0EaJApIBLSikwCDQiKYQjIPAiaIBAQAgmETPv9kbJfRIGkSOqkkv1ZqxZVt27V/d0sVu2655w6hyoLZpSa4nE5VkBMmdG+fnv+3OPPfLjzQ2Zum+k6TsklwpGYUbzeOYcBm2vRILWr60R+YdG9c+l5ch477p5I3V7tXMfxCSsgpkwZ22Us3Zp2Y8ynYzh8+rDrOCVSWmYaQ9c+TqNydYlZepItt4/jwgXXqUq2hPnfEDH/d+yrfTPtZz3lOo7POCkgIlJTRJaLyD7Pvz+bE0FEeojI1ny3DBEZ4HnufRH5Jt9z7X1/FsYfBQYEMmPgDFSVoXFDbdr3X/DHpX/k4KmDzBwynxO9H+P+E68zOzredawS69SxTAKHRBMQAPU/m11qZtotCFdXIOOBz1S1OfCZ5/FPqOoqVW2vqu2BCOAcsCzfLk/++LyqbvVJalMqhFQPYUrfKaw5vIYX17/oOk6J8vHej3ln8zs82eVJujfrzq/m/Z1va3egf9xwPotJch2vxFGFNd3+RIfMjRybNJWqbUNdR/IpVwUkEojx3I8BBlxh/yjgU1U9V6ypTJkxtN1Q7rnhHp5e+TSff/u56zglQkp6CiPjR9K2Xlue7/F83sby5am7MpYKAZlUfOg+Du0vgoXUS5Elj35C5L6X2NL5YZpPiHIdx+dcFZB6qvrjYg1HgXpX2D8amHPRthdEZJuIvCoilxzuICKjRCRBRBJSUlKuIrIpTUSEd/q9Q0j1EKIXRHPy3EnXkZxSVUbGj+R0xmlmDpz5kx8MVmjTnLTJb9Iley3ruo63/hCPbZ8mc+MbwzhQpS3tVrziOo4TxVZARGSFiHz9C7fI/Ptp3njKS46pFJEGQBsg/y/AJgAtgRuBmsC4S71eVd9R1XBVDa9Tx9ZtNv+vWvlqzLt7HsfTjzMsbhi5mus6kjOvfvEq/977bybfNpk29dr87Pl6fxzCwTvGMOTYy8zqd/F3ubLnRFIGOQPuoqJkUGNpLAGVKriO5ESxFRBVvU1VW//CbRFwzFMYfiwQxy/zVoOBhaqale+9v9c8F4D3gE7FdR6mdAtrEMYrt7/CJ/s+4eX1L7uO48SXSV8ybsU4BrYcyKM3PXrJ/X4V9woHG3cjevlIFv+l7HY75mQrX900mg6ZG/n+bzOo2aWl60jOuGrCigeGee4PAxZdZt97uaj5Kl/xEfL6T74uhoymjHjkxkeIahXFhM8msObwGtdxfOqH8z9wz4f30LhaY6b1n3b5adqDg2myYT5p5WrS6umBbF5aNpuEP+77T/oemc6Wfs9w3VN3uY7jlKsC8jegl4jsA27zPEZEwkVk6o87iUgI0AT4z0WvnyUi24HtQG1gkg8ym1JKRJjabyrX1byOqHlRfJv6retIPqGqDF80nCNnjxAbFVugFQaDG9cjOP4j6nOUnH6RJO0774OkJceScau4c8UfSGzWjw5xE13Hcc5JAVHVk6raU1Wbe5q6fvBsT1DVh/Ltd0hVG6n+tHFaVSNUtY2nSWyIqqb5+hxM6XJNhWtYFL2ICzkXGDB3AOeySv+AvxfXv0j8nngm95pMp0YFbwWu0bsTx1+eScesL9jd6QHSz5aNvqON7+3g5skDSarcglabZpbKBaIKy/4Cxni0qN2C2XfNZuvRrTwU/1Cpni9r8b7FjF8xnsE3DOaxmx4r9OubPj6I3SNf4rbTC1gWNo7sUj6698CaZBqO7EtmUCVqrF9McK1qriOVCFZAjMnnzl/fyaSIScz5ek6p/ZHhnhN7uG/BfbSr347p/ad7vTxtq3cfZ1u30Qzc/xIf3fLSZddT92dH954hs9cdVOcUWXGLuaZtM9eRSgwrIMZcZELXCQy+YTDjV4xn/o75ruMUqdSMVCLnRhIcGEzcPXFXN8uuCG1Xvc72VoMZvPFJ4vv8q+iClhAnv03nuw79uC5zJ0mvLaDRnTZrUn5WQIy5iIjwfuT7dG7SmQcWPsDaw2tdRyoSmTmZ3DXvLg6cOsCHd39Is+pF8E06MJDWW2aS2LQfkctGsyT6/at/zxLi7PHzHGzdn7Bz69jz9ExaPnq760gljhUQY35BxeCKxEfH06x6MyLnRrIrZZfrSFclV3MZvmg4K79ZybT+0/hNyG+K7L2lXDCtd84jsV4vesWO5NPomCu/qIRLO5HBrpYD6Hh2FYmPx9D6z/e4jlQiWQEx5hJqVarFkvuXEBwYTO+Zvf16+vcJKyYwe/tsXoh4gaHthhb5+wdWrkCrPXHsqBdB39gHWRH5jyI/hq+cTkpj53X96XRqGRtHTSPslSGuI5VYVkCMuYzQGqEsuX8JZzPPEjEjguQzya4jFdrL619m8vrJPBz+MBO6Tii24wRfU4lWBz5mY+OB3Bb/GKsjnkdz/atn/cSuFJJaRBCWupKE0e9x89vDXUcq0ayAGHMFHRp0YOmQpaSkpxAxI4KjaUddRyqw1754jbHLx3J3q7uZ0neK1yOuCiqocnnC9s1j7bXDuHXVs6xrNYrsc5nFesyi8t26w5xp341rz20n8dmFhL/xoOtIJZ4VEGMKoFOjTnx6/6ckn0kmIsY/rkSmfDmFx5c+zqDrBzHrrlkEBgT65LhBFYLoumc6q2/5E932TGV349tI+6ZkT3uy/Z9rqND9RmpmHWP3P5bTcWI/15H8ghUQYwrolqa3sPj+xSSdSaLL9C7sObHHdaRLenXDqzy65FEGthzInEFzCA4M9unxJTCAW9dN4rOH5nDtqa9IbXEjhz/a5NMMBaLKxgf/RcsxPTkbVJNTn2ygw+9tDfiCsgJiTCF0b9ad1Q+u5nzWebq+15WEIwmuI/2EqjJu+TieWPYEg64fxNyouT4vHvn1fDeaxClrISeHBoM6kzj0ZcgtGVOfnE/+gc3X3k2nmNEk1OpNjb1fEtq37M6s6w0rIMYUUliDMD4f8TlVylXh1vdvZcHOBa4jAZCRncGwuGFMXj+ZR8IfITYqlnKB5VzH4uYx4bA1kfU1f0u7D8ayu2kv0rcdcJrpm7eXkRrSjtbfxLP4N38nPDmeGiHXOM3kj6yAGOOF5rWas37EetrUa0PU/CieWfmM0wWpvkv9ju7vdeeDbR8wqcck3rjjDZ/1eRREozY16XJkAfN7vU2D5AQC27dm34i/4uvlDTMOHWXbDfcS+rvenNXKbJqygTtWP0Vwefso9Iqqlplbx44d1ZiilJGVoSPiRigT0Z4xPTUpNcnnGZbsW6J1JtfRqn+pqnG74nx+/MLaGJesS6oMUgU9VqmZHp0co5qdXazHzD2dqjsHP6tnA6rqecrrgnbPaUpSRrEeszQBEvQXPlOdf6j78mYFxBSH3NxcfXfTu1rphUpa8+81NfbrWM3NzS3246ZdSNOHP35YmYi2+mcr3ZWyq9iPWVTOn1f9YOgy3RwQpgr6ffUWevz5N1XT04v0ODnfJeu+e5/RU0G1VEGXVI3StdP3FukxyoISVUCAu4EdQC4Qfpn9+gB7gP3A+HzbQ4EvPdtjgXIFOa4VEFOc9p7Yq53e7aRMRHvN6KW7U3YXy3Fyc3N1wc4FGvJaiMpE0SeWPKHns84Xy7GK25GkHH27Z6wmSEdV0LPBNfRwn//R7OUrvb8qSUvT1HdjdX/7QZpJkOYgurRCP134dIJmZRVt/rLiUgVE8p7zLRG53lM83gbGqurPhrKISCCwF+gFJAFfAfeq6k4RmQd8pKpzReQtIFFV37zSccPDwzUhoWSNmjGlS3ZuNm8lvMXTK5/mXNY5RnQYwfiu4wmpHnLV762qrPt2Hc+ufpZVh1bRum5r3uj7RpHOa+XKkWTl3+PWUXP+W/TNXEQV0jkXXI0fWnenSp9uXNO5FdKyBdStC1Wr5i3mlJMDZ8/CoUPk7D3AsaVbyP7P59T75gvK52ZwjLosr3M/VcaN5o7fX0s59+MJ/JaIbFLV8J9td1FA/ntwkdVcuoB0Biaqam/P4x/nYPgbkALUV9Xsi/e7HCsgxleOpR1j4uqJTNsyDUUZfMNgRrQfQY/QHgRI4Tps0zPTid8Tz5SNU9iQtIHalWrz3K3PMarjKIICgorpDNzIzIRlC9PZ9/piam5aTufM1fyafT/ZJxchJyCY4Nyf/sI9m0C20IGvq3clq28kN43tRtsOgRTzj+/LBH8sIFFAH/UscSsiDwA3AROBL1T1Os/2JsCnqtr6EscYBYwCaNq0acfDh/13Qjzjf5LOJPHi5y8SkxhD6oVUGldrTO9re9MjpAfhDcMJqR5C+aDyP3lNemY6u0/sZkPSBtYcXsMn+z7hXNY5QquHMrbLWB5s/yCVgis5OiPfycmBxETYvOIHUjfuIWfXXoJOnyD4XCoBOZlkl6tMbsXK5DZuQvnrr6Vet19zS+8qNGzoOnnp4/MCIiIrgPq/8NSfVHWRZ5/VFHMByc+uQIwr57POs2jPImJ3xLL60GpOZ5wGQBDqV6lPucByBAYEcvLcSVIvpP73dQ2rNiSyRSTRraPp2rRroa9ejCkKlyogxXb9q6q3XeVbJANN8j1u7Nl2EqguIkGqmp1vuzElVsXgikS3jia6dTQ5uTkkHktkx/EdHDh1gOQzyWTlZpGdm02NCjVoVK0RodVD6dykM02qNSn2CRCN8VZJbkD9CmguIqHkFYho4D5VVRFZBUQBc4FhwCJ3MY0pnMCAQMIahBHWIMx1FGOuipPrYREZKCJJQGfgExFZ6tneUEQWA3iuLsYAS4FdwDxV3eF5i3HAEyKyH6gFTPP1ORhjTFnntBPd16wPxBhjCu9SfSDWI2eMMcYrVkCMMcZ4xQqIMcYYr1gBMcYY4xUrIMYYY7xiBcQYY4xXytQwXhFJAbydDKs2cKII4/iav+cH/z8Hf88P/n8O/p4f3JxDM1Wtc/HGMlVAroaIJPzSOGh/4e/5wf/Pwd/zg/+fg7/nh5J1DtaEZYwxxitWQIwxxnjFCkjBveM6wFXy9/zg/+fg7/nB/8/B3/NDCToH6wMxxhjjFbsCMcYY4xUrIMYYY7xiBaQARKSPiOwRkf0iMt51nsIQkekiclxEvnadxRsi0kREVonIThHZISKPuc5UWCJSQUQ2ikii5xyec53JGyISKCJbRORj11m8ISKHRGS7iGwVEb9b10FEqovIhyKyW0R2iUhn55msD+TyRCQQ2Av0ApLIWynxXlXd6TRYAYlIdyANmFGQdeNLGhFpADRQ1c0iUhXYBAzwl78/gOStSVtZVdNEJBhYBzymql84jlYoIvIEEA5UU9Xfus5TWCJyCAhXVb/8IaGIxABrVXWqiJQDKqnqaZeZ7ArkyjoB+1X1oKpmkreMbqTjTAWmqmuAH1zn8Jaqfq+qmz33z5K3OmUjt6kKR/OkeR4Ge25+9c1NRBoDdwJTXWcpi0TkGqA7ntVXVTXTdfEAKyAF0Qj4Lt/jJPzsA6y0EJEQoAPwpdskhedp/tkKHAeWq6q/ncNrwFNArusgV0GBZSKySURGuQ5TSKFACvCepxlxqohUdh3KCojxCyJSBVgA/EFVz7jOU1iqmqOq7YHGQCcR8ZvmRBH5LXBcVTe5znKVuqpqGNAXGO1p3vUXQUAY8KaqdgDSAef9sVZAriwZaJLvcWPPNuMjnn6DBcAsVf3IdZ6r4Wl2WAX0cZ2lEG4B+nv6EOYCESIy022kwlPVZM+/x4GF5DVP+4skICnfleuH5BUUp6yAXNlXQHMRCfV0XEUD8Y4zlRmeDuhpwC5VfcV1Hm+ISB0Rqe65X5G8ARm73aYqOFWdoKqNVTWEvP//K1V1iONYhSIilT2DMPA0/dwO+M3IRFU9CnwnIi08m3oCzgeSBLkOUNKparaIjAGWAoHAdFXd4ThWgYnIHOBWoLaIJAHPquo0t6kK5RbgAWC7pw8B4H9VdbHDTIXVAIjxjOgLAOapql8OhfVj9YCFed9HCAJmq+oSt5EK7ffALM8X2YPAcMd5bBivMcYY71gTljHGGK9YATHGGOMVKyDGGGO8YgXEGGOMV6yAGGOM8YoVEGOMMV6xAmKMMcYrVkCMcUhEbhSRbZ41Qyp71gvxm3myTNlmPyQ0xjERmQRUACqSN9/RXx1HMqZArIAY45hnaoqvgAygi6rmOI5kTIFYE5Yx7tUCqgBVybsSMcYv2BWIMY6JSDx506SHkrd87xjHkYwpEJuN1xiHRGQokKWqsz2z9a4XkQhVXek6mzFXYlcgxhhjvGJ9IMYYY7xiBcQYY4xXrIAYY4zxihUQY4wxXrECYowxxitWQIwxxnjFCogxxhiv/B+RauzN53jiCAAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        ""
      ],
      "metadata": {
        "id": "shpAT3p4r4uE"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}
